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[57] ABSTRACT 

An optoelectronic tomographic reconstruction system uti- 
lizes spacial light modulators in charge couple devices to 
perform projection iterative reconstruction techniques and 
simultaneous iterative reconstruction techniques. A back 
projection processor uses a linear array of analog spacial 
light modulators in a cylindrical lens, an image detection 
array, and an image rotator. The image rotator smears a 
projected image at the same angle as the projection was 
taken. The back projection processor thus smears the pro- 
jection back to the image space. An optoelectronic forward 
projection processor uses an spacial light modulator array, 
an image rotator, and an image detecting array. A recon- 
structed image displayed by the spacial light modulator is 
smeared by the rotator to forward project the reconstructed 
image on the image detecting array at the same angle as 
when the measured projection was taken. The forward 
projection processor thus smears the reconstructed image 
back into the projection space. 
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OPTOELECTRONIC SYSTEM FOR 
IMPLEMENTATION OF ITERATIVE 
COMPUTER TOMOGRAPHY ALGORITHMS 

This application is a continuation of U.S. patent appli- 
cation Ser. No. 07/880 ,716, filed on May 8, 1992, now U.S. 
Pat. No. 5,414,623, issued on May 9, 1995. 

TECHNICAL FIELD OF THE INVENTION 

The present invention pertains generally to the field of 
computer tomography, and more particularly to methods and 
apparatus for tomographic reconstruction. 

BACKGROUND OF THE INVENTION 

Computer Tomography (CT) systems have enabled phy- 
sicians to perform noninvasive examinations and glean 
information about the patient by obtaining a cross- sectional 
view of the patient. The quality and reconstruction speed 
have improved dramatically over the years. Simultaneously, 
wider recognition of the value and potential of CT systems 
as well as declining costs have contributed to dramatic 
growth in the industry. Although CT systems have witnessed 
significant increase in market penetration, the untapped 
potential is enormous. In order to tap the market further, both 
costs and performance characteristics have to improve. 
Recent attempts to meet these needs have led, for example, 
to the development of videographic tomographic techniques 
based on the convolution back projection algorithm. 
Although the design of systems that can generate images at 
frame rates (30 frames per second) using the technique are 
now feasible, the method suffers from several disadvantages. 
The disadvantages include: 

High costs due to the need for using expensive high speed 
filters. 

Distortion introduced by the high speed filter, optical 
devices and prism image rotator. 

Commonly used tomographic algorithms can be broadly 
classified as direct and iterative methods, based on the 
method of computation. When a system is underdeterrnined, 
commonly used iterative algorithms, such as the Simulta- 
neous Iterative Reconstruction Technique (SIRT) do not 
guarantee a unique solution. The Projection Iterative Recon- 
struction Technique (PIKT) proposed in the present inven- 
tion is a basic iterative image reconstruction algorithm 
which can be considered as a counterpart of the conventional 
algorithm — SIRT. The PIRT attempts to obtain the 
minimum-norm solution of an under determined system, 
whereas, conventional methods are usually based on an 
attempt to estimate the least squares solution of an overde- 
termined system. Therefore, when the PIRT algorithm is 
applied to an underdetermined system, a unique solution is 
intrinsically guaranteed. The proposed algorithms include a 
family of accelerated algorithms, such as PIRT-CG (PIRT - 
Conjugate Gradient), and PIRT-PC (PIRT - Partial 
Convolution). 

A 2-D cross-sectional image of a 3-D object is shown in 
FIG. 6. Parallel x-ray projections are taken around the object 
in several orientations. The objective of the tomographic 
image reconstruction algorithm is to reconstruct the 2-D 
cross-sectional image on the basis of information contained 
in the ray- sums of projections measured from several ori- 
entations across the image plane around the object. 

In the continuous case, a ray sum is expressed by the 
Radon transform which is obtained by performing a line 
integration along each ray. As shown in FIG. 7, a line 
integral P e (t) can be defined as 



Pe(t) 



■ f > 

J (9, /) line 



(1) 



5 Using the sifting property of the delta function, eq. 5 can be 
rewritten as 



10 



(2) 



yfc;y)5(*cose + ysinG - t)dxdy 



The function P e (t) is known as the Radon transform of the 
function f(x, y). 

In the discrete case, eq. 6 can be rewritten in the form of 
15 a summation and projection operations can be modeled 
using a linear system representation 



Ax=b 



(3) 



where x represents all pixels on the two dimensional image, 

20 b represents data measured at all projection orientations, and 
the matrix A maps data from the image space to the 
projection space. Image reconstruction involves estimation 
of the 2-D image x from known projections b. 

Tomographic image reconstruction involves operations of 

25 mapping the data from the projection space back to the 
image plane. The operation is called back projection. The 
back projection was one of the earliest methods used to 
obtain the a cross-section of an object in an x-ray film before 
computed tomography was invented. Since the method 

30 simply smears the projection data back to the image space 
instead of solving for the true inverse, information in the 
reconstructed image severely lacks in detail. 

Tomographic image reconstruction involves estimation of 
f(x, y) from given P e (t) using the inverse Radon transform. 

35 The task is summarized in eq. 6 for the continuous case. In 
the discrete case, the objective is to determine x from known 
b as expressed in eq. 7. In computed tomography, the 
enormous amount of data contained in the projections col- 
lected in several directions has to be appropriately manipu- 

40 lated to obtain the spatial distribution of the parameters. The 
algorithms for tomographic imaging solves the inverse prob- 
lem by estimating the cross-sectional images from the given 
projections. 

In the case of most iterative reconstruction methods, the 

45 error corrections are fed back in the image space except for 
a few exceptions, such as the Projection Space Iterative 
Reconstruction-Reprojection (PIRR) and Projection Space 
MAP (PSMAP) methods. The PIRR projects the recon- 
structed image to the projection space recursively in order to 

50 recover the missing projection data whereas the PSMAP 
optimizes data in the projection space iteratively and then 
reconstructs the image using convolution back projection 
(CBP). The approaches and objectives of the projection 
space iterative algorithms differ from the PIRT proposed in 

55 the present invention. 

Because of the extreme computational demands, iterative 
methods usually are not able to compete with direct algo- 
rithms in commercial CT systems. However, the iterative 
algorithms still offer advantages in certain applications. 

60 They are particularly suitable for reconstructing images 
from incomplete data, reconstruction with a priori statistical 
knowledge as well as single photon emission computerized 
tomography (SPECT) and positron emission tomography 
(PKT). 

65 As mentioned before, conventional x-ray tomography 
techniques were used to estimate cross- sectional images of 
objects even before the invention of computed tomography. 
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In this method, a photographic film and the object are rotated 
synchronously. X-rays pass through a narrow slit, penetrate 
the object and are recorded on the film. The signal is then 
smeared on the film plane. Although, the technique may be 
considered as the earliest optical implementation of 5 
tomography, the method is mathematically equivalent to the 
back projection method. However, the quality of results 
obtained are poor compared with those obtained using 
modern computed tomography techniques. 

The computed tomography methods differ from the con- 
ventional tomography schemes in that they attempt to solve 
the inverse problem using the measured projection instead of 
simply smearing the projections back into the image space. 

The earliest reported optical computed tomographic 
reconstruction processor was built by Peter in 1973. An 
image was first recorded on a film using back projections. 15 
The output image was obtained by filtering the blurred 
image using a coherent optical spatial filter. The resulting 
image was much sharper than those obtained using back 
projections only. 

In order to avoid the problems associated with coherent 20 
processing such as speckle and other coherent noise, several 
incoherent optical tomographic reconstruction systems were 
built. 

The Oldelft transaxial tomography system, built in 1978, 
implements the convolution back projection method. The 25 
Oldeft system is a hybrid system and optics is used only for 
one dimensional and two channel convolutions. The two 
channels are used for positive and negative valued convo- 
lution respectively. 

The Edholm's system, built in 1977, is an optical system 30 
using films. The original projections and the filtered negative 
projections are prepared on two separate films. The recon- 
structed image is then recorded on a rotating output film. 

Since 1977, several structures have been proposed by 
Gmitro et al. where pupil plane masks have been used for 35 
spatial radius filtering operations. This approach, known as 
optical transfer function (OTF) synthesis, is a technique for 
performing spatial filtering operations in an incoherent sys- 
tem. The loop processor records all projections on a con- 
tinuous film loop and the drum processor records projections 40 
on the surface of a drum The CCD processor collects the 
back projected output image using a CCD camera. 

Several coherent computed tomography systems have 
been proposed by Hansen et al., and Nishimura and 
Casasent These approaches have been summarized by 45 
Gmitro et al. All of these algorithms involve the use of direct 
algorithms. Due to the finite dynamic range of materials and 
devices, and distortions of optical transforms, these 
approaches could not compete with electronic computers in 
respect of the quality of reconstructed images. 50 

Advances in technology related to video imaging devices 
have led to improvements in the quality and speed of optical 
implementations. Recently, a videographic tomographic 
structure for medical imaging, built by Gmitro et al., was 
able to achieve 1% contrast resolution. The structure was 55 
able to achieve real time reconstruction. Three filter struc- 
tures were evaluated including a digital FIR filter, an 
Acoustic-optic (AO) convolver and a Surface Acoustic 
Wave (SAW) convolver. The best results were obtained 
using the digital FIR filter. However, the digital FIR filter 60 
was not only expensive but also introduced distortions in the 
low frequency range since the order of the filter used was not 
long enough to cover the entire frequency range. In addition, 
optical distortions were not eliminated. Nevertheless, the 
development is very encouraging since it demonstrated the 65 
feasibility of high quality and high speed optoelectronic 
tomography. 
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The present invention provides a new approach/structure 
for tomographic reconstruction. The basic structure is built 
from optoelectronic devices and allows for implementation 
of a number of reconstruction algorithms. The critical 
devices used in this structure are Spatial Light Modulator 
(SLM) and Charge Coupled Device (CCD) arrays which are 
commonly used as liquid crystal television display panels 
and solid state video image detectors, respectively. 

Spatial Light Modulators (SLM's) are devices which can 
be employed to modulate the intensity, magnitude, polar- 
ization or phase of light. Applications of SLM range from 
commercial television displays, real-time image processing 
to parallel optical computing. FIG. 8 shows an example of 
a linear SLM array with four cells. The optical transmissiv- 
ity of each cell can be controlled by the applied modulating 
signal. The intensity of the output light beam from each cell 
is, therefore, a function of the intensity of the incident beam 
and the transmissivity of the cell. SLM's can be classified on 
the basis of their addressing modes. If the modulating signal 
is controlled by an electrical signal, it is referred to as an 
electrically addressed SLM (E-SLM). If the modulating 
signal is a second beam of light, it is classified as an optically 
addressed SLM (O-SLM). 

SLM's have been built using several technologies. This 
had led to the development of the optoelectronic SLM, 
opto-acoustic SLM, and opto-magnetic SLM. Commercially 
available SLM's, such as Liquid Crystal (LC) SLM, and 
Ferroelectric Liquid Crystal (FLC) SLM will be briefly 
described in this section. For the sake of completeness, 
Multiple Quantum Well (MQW) SLMs are also discussed. 

Liquid crystals are organic materials that possess an 
intermediate phase between the solid and liquid phase. The 
molecular orientations of liquid crystal materials can be 
changed by applying electrical fields. Therefore, the polar- 
ization of the light through such materials can also be 
twisted. Amplitude and intensity modulation is obtained by 
placing a polarizer and an analyzer in front of and behind the 
liquid crystal layer respectively. The technology relating to 
Liquid Crystal SLM has been used in commercial video 
image projectors and miniature television sets. Since com- 
mercial LCTV's offer many of the same attractive features 
as other modulators, but at only a fraction of the cost, they 
have also been used in many optical signal processing and 
computing systems. 

Ferroelectric liquid crystals are characterized by a spon- 
taneous molecular polarization caused by an anisotropy in 
the molecule. The molecular polarization allows the orien- 
tation of FLC's to be easily switched with a small electric 
field. The FLC SLM has the advantage of high speed and 
high contrast ratio. Commercial available FLC SLM only 
offer binary light modulation since tilting of FLC molecules 
is confined to two orientational positions. 

Multiple quantum well (MQW) structures consist of thin 
layers of low bandgap semiconductor (wells) sandwiched 
between layers of larger bandgap semiconductor (barriers). 
When the thickness of the well layers is on the order of a 
carrier de Broglie wavelength, the electron and the hole are 
forced to orbit close to each other and the binding energy 
increase correspondingly. The electrical and optical proper- 
ties of the structure are then dominated by quantum size 
effects (QSE). The QSE results in the features of step-like 
absorption edges in the optical absorption spectrum and 
room temperature exciton resonances. The change in the 
energy levels of the excitons resulting from applied electric 
field, called the quantum-confined Stark-effect (QCSE), 
allows shifting of the abrupt, highly absorbing edge. By 
shifting the absorption edges, the MQW structures produce 
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larger absorption changes (in a narrow spectrum range 
around the absorption edges) than those in bulk semicon- 
ductors with the same applied field. 

The LCTV's are able to provide better grey level images. 
However most of them can only operate at television frame 5 
rates (30 to 60 frames per second). The FLC SLM's support 
binary processing only but can operate at relatively higher 
speed (up to 100 Khz). The MQW SLM's are expected to 
operate at GHz rates. However, the spectrum of modulated 
light has to be within a narrow range. 



TABLE 1 



Characteristics of Several Common SLMs 



SLM 


Visibility 


Resolution 


Size 


Speed 


Cost 


MSLM 


0.5 


4 lp/mm 


25 mm dia 


2 sec 


$25K 


MOD 


0.91 


6.4 lp/mm 


1 x 1 cm 


200 Hz 


$18K 


DMD 


0.5 


10 lp/mm 


0.64 cm 


60 Hz 


? 


FELC 


0.9 


40 lp/mm 


12.5 Tnm dia 


60 Hz 


$17.5K 


HCLCLV 


0.86 


60 lp/mm 


50 X 50 mm 


60 Hz 


$25K 


Epson 


0.96 


6.3 lp/mm 


2.54 x 1.9 cm 


60 Hz 


$800 


LCTV 













TABLE 2 



Specificaitons of Several Electrically Addressed SLMs 



frame 
rate 



pixel 



fill contrast 



device 


material 


pixels 


Hz 


size pm 


factor 


ratio 


STC Ltd. 


FLC 


128 x 128 


165 


165 


x 165 


0.83 


200:1 


Displaytech 


FLC 


10 x 10 


2000 


l 3 


x 10 3 


0.77 


100:1 


Disp.tech 


FLC 


64 x 64 


4500 


45 


X45 


0.56 


12:1 


CMOS 
















Semetex SMD 


Iron 


128 x 128 


2000 


56 


x56 


0.54 




Litton 


Iron 


128 x 128 


100 


56 


x56 


0.54 




MOSLM 


Garnet 














TI DMD 


Def\mirr 


128 x 128 


1200 


25 


x 25 


0.9 


2:1 




or 
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Charge Coupled Devices (CCD) are also called Charge 40 
Transfer Devices (CTD). These devices were introduced by 
Boyle and Smith in 1970. The applications of CCD's include 
optoelectronic computing, charge domain signal processing, 
focal plane image processing, high speed analog-to-digital 
and digital-to-analog conversion, time-axis conversion, and 
image detection. 

A CCD array functions like shift registers in which 
sampled values of an analog signal are stored in the form of 
charges in a series of neighboring cells. Clock pluses allow 
the transfer of charge from one cell to the next without 
significant loss in accuracy. The charge, which is called a 50 
charge packet, is a small amount of charge stored in potential 
wells created by the gate voltage. By periodically varying 
the electrode or gate voltage, the potential wells are shifted 
along the semiconductor. The charge packet, located under 
the voltage gate, moves along with the potential wells. 55 

CCDs can be classified as surface-channel CCD (SCCD) 
and bulk-channel CCD (BCCD) depending on the device 
structures. A device is called SCCD if the charge resides at 
the semiconductor surface. In the BCCD, the charge location 
is moved away from the surface into the n-channel. There is 60 
no interaction between the charge and the interface states. In 
the case of peristaltic CCD, the n-channel is thicker and the 
charge is farther away from the surface. In spite of the 
increased process complexity, almost all of today's devices 
are BCCD's, due to their superior performance. 65 

High performance Silicon CCDs can offer charge transfer 
efficiencies as high as 0.999999. Some of the fastest Silicon 



CCDs have been operated at several hundred MHz. Silicon 
peristaltic CCDs have been operated with clock frequencies 
up to 200 MHz. High mobility GaAs is a good material for 
building very high frequency CCD's. Devices with operat- 
ing frequency in the range of GHz have been fabricated. 
However, typical operational frequencies of commercially 
available Silicon CCD's are usually below 20 MHz. 

CCD detecting arrays are commonly used as solid state 
image detectors. FIG. 9 shows a linear CCD detector array 
with four cells. Each cell generates an electrical charge 
proportional to the number of photons incident on it. 
Alternatively, the charge developed is proportional to the 
intensity of the light integrated over the period of exposure. 
Charges in a CCD array can be shifted from cell-to-cell 
without significant loss in magnitude. The output is a 
sequence of voltage values proportional to the charge in each 
cell. Typical contrast resolution of commercial image detec- 
tors can be in excess of 4096 distinguishable gray levels. 

SUMMARY OF THE INVENTION 

The present invention provides for iterative tomographic 
reconstruction. The apparatus of the invention comprises an 
optoelectronic back projection processor, an optoelectronic 
forward projection processor, and a projection space auxil- 
iary processor. The system optionally includes an image 
space auxiliary processor. 

The back projection processor comprises a linear array of 
analog SLMs, an image detection array, and an image 
rotator. A light source is applied to the analog SIM array, 
which transmits an image from the projection space through 
the image rotator back to the image detecting array. The 
image rotator smears the projection image, initially obtained 
by a parallel beam imaging device, on the image detecting 
array, at the same angle as the projection was taken. The 
back projection processor thus smears the projection back to 
the image space. 

The optoelectronic forward projection processor com- 
prises an SLM array, an image rotator, and an image 
detecting array. A light source passes through the SLM array 
and the image rotator and is projected on the image detecting 
array. A reconstructed image displayed by the SLM is 
smeared by the rotator to forward project the reconstructed 
image on the image detecting array at the same angle as 
when the measured projection was taken. The purpose of the 
forward projection processor is to smear the reconstructed 
image back into the projection space. 

The projection space auxiliary processor receives data 
from the forward projection detecting array, performs simple 
arithmetic operations, stores projections, and sends data 
back to the back projection SLM array. The projection space 
auxiliary processor is used to perform successive error 
corrections using back and forward projections. 

According to another aspect of the invention, a projection 
iterative reconstruction technique (PIRT) uses the above- 
described structure to reconstruct an image by reiteratively 
correcting a plurality of state vectors representing the state 
of the reconstructed image. The state vectors for each 
projection angle are supplied to the SLM array in the back 
projection processor, which uses the state vectors to recon- 
struct an iteration of the image on the CCD detecting array. 

The approach of the present invention employs feedback 
to overcome distortion introduced by the optical devices. 
The system is based on existing technology and offers 
significant improvement over current systems both with 
respect to cost and performance. The advantages offered by 
this invention are: 
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a) It is able to compensate for distortions caused by optical 
systems. 

b) It is able to reduce effects caused by the limited dynamic 
range of optoelectronic devices. 

c) It does not require computation of Fourier transforms or 
convolution operations. 

d) It is a closed-loop system. 

e) Since all computations are performed in the projection 
space, the effects of dark current in the CCD detector do 
not result in degradation of the reconstructed image. 
Since there is no need for computing Fourier transforms 

and convolution operations, this structure can be used for 
high speed and low cost tomography reconstruction where 
the rotations can be replaced by static geometric rotation 
schemes. Static image rotation techniques can be performed 
using devices employing current technologies, such as holo- 
grams or coupled CCD/SLM pairs. The scheme described in 
the disclosure, however, employs mechanical image rota- 
tors. The operation of the system is independent of the 
technique employed for rotating the image. 

Iterative tomographic reconstruction algorithms that can 
be implemented/executed using the structure include: 

a) Projection Iterative Reconstruction Technique (PIRT). 

b) accelerated version of PIRT using partial convolutions 
(PIKT-PC). 

c) accelerated version of PIRT using conjugate gradient 
acceleration technique (PIRT-CG). 

d) Simultaneous Iterative Reconstruction Technique (SIRT). 

e) accelerated version of SIRT using filtered back projection 
method (SIRT-FBP). 

BRIEF DESCRIPTION OF THE DRAWING 

FIG. 1 illustrates a block diagram of the Projection 
Iterative Projection Technique (PIRT) according to the 
present invention; 

FIG. 2 is a block diagram of the Simultaneous Iterative 
Reconstruction Technique (SIRT); 

FIG. 3 is a block diagram of the optoelectronic architec- 
ture for iterative image reconstruction according to the 
present invention; 

FIG. 4-a is an illustration of two projections of a 3-D 
object; 

FIG. 4-& is an illustration of the back projection processor 
according to the present invention; 

FIG. 5 is an illustration of the forward projection proces- 
sor according to the present invention; 

FIG. 6 depicts a two dimensional cross- sectional image of 
a three dimensional object; 

FIG. 7 depicts a two dimensional cross-sectional image 
and corresponding projections; 

FIG. 8 depicts a linear spatial light modulator (SLM) 
array with four cells; 

FIG. 9 depicts a linear charge coupled device (CCD) 
detecting array with four cells; 

FIG. 10 depicts a two dimensional cross-sectional image 
and the corresponding projections; 

FIG. 11 depicts a least squares solution for an overdeter- 
mined system; 

FIG. 12 depicts a block diagram of a simultaneous itera- 
tive reconstruction technique (SIRT); 

FIG. 13 depicts an optoelectronic forward projection 
processor using a one dimensional SLM array and a two 
dimensional CCD array; 

FIG. 14 depicts an optoelectronic back projection proces- 
sor using a two dimensional SLM array and a one dimen- 
sional CCD array; 
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FIG. 15 depicts an iterative optoelectronic structure 
implementing SIRT; 

FIG. 16 depicts a nonlinear response curve of a GaAs/ 
AlgaAs multiple quantum well SLM; 
5 FIG. 17 depicts a block diagram for an optoelectronic 
implementation of SIRT with nonlinear functions of forward 
projection and back projection SLM arrays; 

FIG. 18 depicts a Shepp and Logan phantom consisting of 
ten ellipses; 

FIG. 19 depicts a cross- section of the original phantom 
and actual reconstructed images; 

FIG. 20 depicts a cross-section of displayed images after 
contrast stretching; 
15 FIGS. 21a-d depict a reconstructed image obtained using 
the SIRT at the 16 th , 32 nd . 64* and 256** iterations, respec- 
tively; 

FIGS. 22a-d depict reconstructed images at 16 th , 32 nd , 
64 th , and 128 th iterations, respectively where the number of 
20 distinguishable grey levels of the forward projection SLM is 
set at 256; 

FIGS. 23a-d depict reconstructed images at the 16* 32™*, 
64* and 128* iterations, respectively where the number of 
distinguishable grey levels of the forward projection SLM is 
25 set at 1024; 

FIG. 24 depicts a plot of residual versus the number of 
iterations for SIRT; 

FIGS. 2Sa-d depict reconstructed images at the 16*, 32 nd , 
30 64 th , and 5 12 th iterations, respectively, with additive noise in 
the image space; 

FIG. 26 depicts a minimum-norm solution for an under- 
determined system; 

FIG. 27 depicts a block diagram of a projection iterative 
35 reconstruction technique (PIRT); 

FIG. 28 depicts projection A and B covering the entire 
region of solution X; 

FIG. 29 depicts a convex support region constrained by 
support regions of projections; 
40 FIG. 30 depicts a conjugate gradient method converging 
in N steps; 

FIG. 31 depicts a steepest descend method converging in 
an infinite number of steps; 
45 FIG. 32 depicts an optoelectronic structure for the PIRT; 
FIGS. 33a-d depict reconstructed images using PIRT at 
the 8* 16* 32"* and 64 th iterations, respectively; 

FIG. 34 depicts a plot of residuals of reconstructed images 
using the original phantom using PIRT and CG-PIRT; 
50 FIGS. 3Sa-d depict reconstructed images using the PIRT- 
CG at the 8*, 16* 32 nd , and 64 th iterations, respectively; 

FIGS. 36a-d depict reconstructed images using PIRT 
with additive noise at the 16*, 32"* 64 th and 512* 
iterations, respectively; 
55 FIG. 37a depicts a projection taken on a x-y plane of a 
cross-sectional image by line integration in a direction 
perpendicular to the x' axis and FIG. 37b depicts the Fourier 
transform of the projection of FIG. 37a; 

FIG. 38 depicts a spatial frequency response of a radius 
60 filter; 

FIG. 39 depicts a block diagram of an iterative filtered 
back projection (IFBP); 

FIG. 40 depicts a two lens system implementing a two 
65 dimensional Fourier transform spatial filter; 

FIG. 41 depicts a frequency response of the coherent 
transfer function for a circular aperture; 
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FIG. 42 depicts a frequency response of the optical 
transfer function for a circular aperture; 

FIGS. 43a-b depict the frequency responses of negative 
and positive channels of a complimentary filter; 

FIGS. 44a-c depict characteristic curves of a negative 
response SLM and a positive response SLM; 

FIG. 45 depicts structure of an optoelectronic implemen- 
tation of an iterative filtered back projection technique; 

FIGS. 46a-d depict reconstructed images obtained using 
the IFBP at the 1* 2 nd , 3 rd and 64 th iterations, respectively 
with finite numbers of distinguishable grey levels; 

FIGS. 47a— d depict an original 32x32 phantom and the 
reconstructed images using the IFBP at the 1**, 256 th and 
30,720 r * iterations, respectively; 

FIG. 48 depicts an impulse response of a radius filter; 

FIG. 49 depicts a line diagram of a proposed hybrid 
prototype; 

FIG. 50 depicts projection data back projected into an 
image space at the same angle as when an original measured 
projection was taken; 

FIGS. 51a-c depict a dove prism rotated 0, 45 and 90 
degrees counterclockwise, respectively; 

FIG. 52 depicts a forward projection processor wherein 
the reconstructed image is forward projected into projection 
space at the same angle which the original measured pro- 
jection was taken; 

FIGS. 53a-b depict a timing diagrams for the control of 
exposure period for an analog SLM; 

FIG. 54 depicts a circuit implementation of an analog 
SLM array using a binary SLM array; 

FIG. 55 depicts four SLM cells and 16 grey levels; 

FIGS. 56a-d depict reconstructive images at the S th y 16 th , 
32 nd and 64 th iterations, respectively obtained using PIRT- 
PC; 

FIGS. 57a-d depict reconstructed images at the S th , 16 th , 
32 nd and 64 th iterations, respectively obtained using basic 
PIRT; 

FIG. 58 depicts cross-sections of reconstructed images at 
the 1** 2 nd , 4 th and S ,h iterations obtained using PIRT-PC; 

FIG. 59 depicts cross-sections of reconstructed images at 
the 1 st , 2 nd , 4 th and % tH iterations obtained using basic PIRT; 

FIG. 60 depicts cross-sections of reconstructed images at 
the 32 nrf iteration obtained using PIRT-BC with constraints 
and assuming finite dynamic range; 

FIG. 61 depicts cross-sections of reconstructed images at 
the 64 th iteration obtained using PIRT-PC with constraints 
and assuming finite dynamic range; 

FIGS. 62a-** depict reconstructed images at the S ,h , 16 th , 
32 nd and 64 th iterations, respectively, obtained using PIRT- 
PC with constraints and assuming finite dynamic range; 

FIG. 63 depicts a cross-section of an image where all 
projections are biased to non-negative values; and 

FIG. 64 depicts cross-sections of two reconstructed 
images where the images are reconstructed by processing 
positive and negative protection separately. 

DETAILED DESCRIPTION OF THE 
INVENTION 

In the following detailed description of the preferred 
embodiments, reference is made to the accompanying draw- 
ings which form a part hereof, and in which is shown by way 
of illustration, specific embodiments in which the invention 
may be practiced. It is to be understood that other embodi- 
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ments may be utilized and structural changes may be made 
without departing from the scope of the present invention. 

Overview 

5 FIG. 1 and FIG. 2 show feedback block diagrams asso- 
ciated with the basic PIRT and SBRT algorithms respectively, 
to be discussed below. 

Optoelectronic Structure for Tomographic 
10 Reconstruction 

FIG. 3 shows a line diagram outlining the basic building 
blocks in the system. The building blocks include 

An Optoelectronic Back Projection Processor. 
15 An Optoelectronic Forward Projection Processor. 

An Image Space Auxiliary Processor. 

A Projection Space Auxiliary Processor. 

Each of the building blocks is described briefly in the 
2Q following sections. 

The Optoelectronic Back Projection Processor 

A set of back projections in tomographic image recon- 
struction involves mapping the data in the projection space 
25 into the image space from all the measured projection 
angles. FIG. 4-a shows the projection of x-rays onto the 
projection space at two different angles. The optoelectronic 
back projection processor as shown in FIG. 4-b performs the 
tomographic back projection. It consists of: 
30 a linear array of analog SLMs and a cylindrical lens. 

an image detection array. 

an image rotator. 

The transmissivity of the linear SLM array represents the 
magnitudes of the data in the projection space corresponding 

35 to a projection angle. The pattern in the linear array is then 
"stretched" into an array of strips as shown in FIG. 4-b. The 
linear SLM array is not specified. However, based on 
available devices, an implementation using binary Ferro- 
electric Liquid Crystal (FLC) SLMs is suggested. 

40 The image detection array integrates the projection data 
corresponding to the area intercepted on the pixel by the 
light beam at all the projection angles. The image detection 
array can be a CCD image detector or a commercial CCD 
TV camera. 

45 A complete set of back projections involves performing 
back projection operations at all the projection angles over 
an angle of 180°. The purpose of the image rotator is to 
rotate the projection pattern to the corresponding projection 
angle relative to the image detector. It can be done by 

50 rotating the linear SLM array and the lens back and forth 
over a range 180°, or by using a prism image rotator driven 
by a stepper motor. Alternative approaches using a static 
structure can also be employed. 

55 The Optoelectronic Forward Projection Processor 

The optoelectronic forward projection processor as shown 
in FIG. 5 performs the tomographic forward projection to 
map the data in the image space into the projection space. 

It consists of: 

an image display array. 

a linear projection detector array and a cylindrical lens for 
"compressing" a two dimensional image to a one 
dimensional image. 
65 an image rotator. 

The image display array displays a reconstructed image 
using a two dimensional analog SLM or commercial liquid 
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crystal projection TV. The linear projection detector array 
detects the image generated by the cylindrical lens. The 
linear projection detector array can be a linear CCD array. 
The role and the design of the image rotator is identical to 
the one described. 5 

The Image Space Auxiliary Processor 

Since operations of image addition can be performed 
directly on the CCD array, a separate image auxiliary 
processor may not be required However, when performance 
of devices is not guaranteed, additional flexibility and 
improvement in accuracy can be gained by employing an 
auxiliary processor. The processor is required to be able to 
store a frame of the image and perform simple pixel by pixel 
addition or subtraction on images. Real-time digital video 
image processors which are commercially available can be 
used to perform these operations. 



The Projection Space Auxiliary Processor 

The basic functions of this auxiliary space processor are: 
to receive data from the forward projection detector array, 
to perform simple arithmetic operations, 
to store projections; and 

to send data to the back projection linear SLM array. 

This basic auxiliary processor can be built using either a 
commercial digital microprocessor or special purpose pro- 
cessor. 

Algorithms 

The rationale behind the algorithms is that a cross- 
sectional image can be reconstructed by successive error 
corrections using back and forward projections. The tomo- 
graphic forward projection operations can be expressed by a 
matrix A and the corresponding back projection operations 
are expressed by the transpose of the matrix, A T . 

PIKT: Projection Iterative Reconstruction Technique 

The present invention provides a new projection iterative 
reconstruction technique (PIKT). The operation involved in 
reconstruction using this technique can be summarized as 
follows: 

fjfc=^i+ot [b—A Xk _ x ] (4) 
A T t k 

where x^: an Nxl vector, represents the reconstructed image 
at the k* A iteration. 



50 



an Mxl vector, represents the measured 
projections. 

an MxN matrix, maps an image into the 
corresponding projections. 
A state vector, is interactively built up in the 
projection space. 



a and p are scalars which are adjusted to achieve proper 
system gain. The block diagram is given in FIG. 1. 

For an underdetermined system, the PIRT is guaranteed to 
converge to the unique minimum-norm solution. 



lim x k =A T [AA T ]~ 1 b 



60 



(5) 
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processor, t K _ x is updated to obtain t k and the result is then 
stored in the projection space auxiliary processor. The 
projection space auxiliary processor is also used to store b. 
The other matrix-vector multiplication operation involved in 
computing p A r t k is performed using the optoelectronic 
back processor to x k which represents the output image. 

PIRT-PC: Projection Iterative Reconstruction 
Technique with Partial Convolution 

When the Convolution Back Projection (CBP) method is 
used, full length convolutions have to be performed. In order 
to reduce hardware cost and accelerate the convergence of 
the PIRT, a short length convolution can be used as an 
acceleration technique. The improved algorithm can be 
described by 



(6) 



20 



or 



(7) 



In eq. 1, the matrix- vector multiplication A x k _ x is per- 
formed using the optoelectronic forward projection 



t k =t^+a H {b-A x k ^] 
x k ^A T t k 

r^r^-a A A T H r k _ x 
**P**-rHX A T Hr k 

where H is an NxN matrix representing a FIR filter of 
reduced length, r* is the residual in the projection space and 
r 0 =b. The order of the filter is significantly less than the order 
of the FIR filter used for the CBP algorithm and typically 
consists of only a few taps. 

Other Techniques 

Applications of the invention involving PIRT-CG, SIRT 
and SIRT-FBP are disclosed hereinbelow. 

Methods for Eliimnating Effects of Limited 
Dynamic Range 

There are three sources which affect the dynamic range of 
reconstructed image: 

a) The non-negative offset associated with the projection 
data may reduce the actual dynamic range to approxi- 
mately 10%. 

b) The finite dynamic range of the back projection 1-D SLM 
array. 

c) The finite dynamic range of the forward projection 2-D 
SLM array. 

The first effect can be eliminated by back projecting the 
positive and negative projection data separately. In addition, 
the error caused by the finite contrast ratio can also be 
canceled out using the dual projection scheme. The iterative 
procedure described in eq. 4 provides a choice for eliminat- 
ing distortions caused by the finite dynamic ranges of the 
SLM arrays. In the iterative procedure, only back projec- 
55 tions and forward projections of residuals are involved. In a 
stable feedback system, the magnitude of the residual will 
vanish as the number of iterations increase. Therefore, the 
vanishing residuals can always be adjusted to fit the maxi- 
mum dynamic range of the SLM arrays and then reseated 
back appropriately using the auxiliary processors after opti- 
cal projections. 

Background on SIRT 
The Simultaneous Iterative Reconstruction 

65 

Technique (SIRT) uses an iterative method to solve linear 
equations for tomographic image reconstruction. The SIRT 
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was first proposed by Gilbert in 1971. The basic SIRT uses 
Richardson's method to solve a system Ax=b. The solution 
x=[A T A]"" 1 A r b represents the least squares solution of the 
normal equations. It can also be treated as a method of 
quadratic optimization. Acceleration techniques, such as the 
Conjugate Gradient method, have been used in conjunction 
with SIRT to expedite convergence [18, 31-32]. 

Even though iterative tomographic algorithms have been 
discussed extensively during the last two decades, most 
applications employ direct reconstruction methods, such as 
the Filtered Back Projection (FBP) method and the Convo- 
lution Back Projection (CBP) method. The computation cost 
of iterative algorithms is formidable if they are implemented 
using today's electronic computers. An optoelectronic struc- 
ture implementing SIRT is presented in this section. The 
proposed structure performs back and forward projection 
operations by simply projecting data optically. Time con- 
suming matrix operations are, therefore, eliminated and the 
processing speed is improved significantly. Since, there is no 
need for computing Fourier Transforms or performing con- 
volution operations explicitly, the cost of such a system is 
significantly lower. 

This section introduces the notation used in casting the 
tomographic reconstruction issue as a problem of solving a 
linear system The SIRT is reviewed as an iterative method 
for solving overdetermined systems. This is followed by a 
discussion on the issues of convergence. The optoelectronic 
implementation of SIRT is introduced after explaining the 
basic building blocks employed for back and forward pro- 
jections. 

A 2-D cross-sectional image of a 3-D object is shown in 
FIG. 10. The objective of the tomographic image recon- 
struction algorithm is to reconstruct the 2-D cross-sectional 
image on the basis of information contained in the ray-sums 
of projections measured from several orientations across the 
image plane around the object. A ray-sum is an integration 
or a summation of pixel grey levels across the image plane 
along that ray. A projection is a set of ray-sums in parallel 
with the same orientation. 

Consider a 2-D discrete model for the nxn cross-sectional 
image where the pixels and ray-sums are numbered from 1 
to N and 1 to M respectively. Let x y , j=l, 2, . . . , N, denote 
the grey level of the ] th pixel, and b„ i=l, 2, . . . , M, represent 
the i th ray-sum Then, the ray-sum b, can be expressed as 



N 

bi — £ WifXj, for i = 1, . . . ,M 
j=l 



where 



b=[b 1 b 2 . . . b M ] r 



In this expression, all the nxn pixels in the 2-D cross- 
sectional image are stacked into the Nxl vector x and all the 
ray-sums from all projection orientations are denoted by the 
Mxl vector b. The matrix A is an MxN projection operator 
matrix mapping x into b and w lV is its (i element. Image 
reconstruction involves estimation of x from known b. 

The operation of mapping the data from the projection 
space back to the image plane is called back projection. The 
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(8) 



where w^- is the weight factor representing the fractional area 
intercepted by the i th ray-sum and the ] th pixel. The projec- 50 
don operations of eq. 8 can be rewritten in a matrix form as 
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back projection operation can also be expressed in the form 
of summation or matrix vector multiplication as 



i=l 



PA 7 *? 



(10) 



(11) 



where p is a constant used to normalize back projections. 
Usually, p is chosen as 



Nli Ml- 



(12) 



where and denote matrix one-norm and infinite- 
norm. 

In what follows, the least squares solution of the normal 
equations representing overdetermined systems is reviewed 
briefly. The Simultaneous Iterative Reconstruction Tech- 
nique (SIRT) for cornputing the least squares solution of 
overdetermined systems is then discussed. 

When M>N and all N columns of the MxN matrix A are 
linearly independent, the system is overdetermined. This 
implies, that there is no direct solution, since the number of 
equations exceeds the number of unknowns. However, there 
exists a unique least squares solution. Using the normal 
equation method, the least squares solution for the overde- 
termined system is 
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^AT 1 A T b 



(13) 



The result can be derived by minimizing ||Ax-b|| 2 where ||.|| 2 
is the Euclidean norm and A T A is symmetric and positive 
definite, 
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(14) 
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-Z-WAx-bW? = -2-iAx-bnAx-b) 

= A r Ax-A T b 
= 0 

A T Ax = A T b 
x = [A T A]- l A T b 
if [A T A\ is nonsingular. 

FIG. 11 shows the least squares solution for an overdeter- 
mined system where M=3, N=2. 

In many applications such as the image reconstruction, 
the dimensions of the matrix A are large. In this case, least 
squares solution of the normal equations cannot be solved 
for directly using matrix inversion. Since the matrix A is 
extremely large and sparse in such cases, the matrix A cannot 
be formed in practice due to memory limitations and con- 
sequently only the non-zero entries of the matrix A are 
calculated on-line on an "as needed" basis and disposed off 
shortly thereafter. Since the matrix is not available in its 
entirety, techniques, such as Singular Value Decomposition 
(SVD), cannot be applied. However, these systems can be 
solved for using iterative methods where the matrix A is not 
required in full. 

In eq. 13, the least squares solution x of the overdeter- 
mined linear system can be obtained by solving for x in a 
linear system D x=c, where D=A T A is an NxN symmetric 
positive definite matrix and c=A r b is an Nxl vector. The 
least squares solution of the overdetermined system can be 
obtained using iterative methods, such as die basic RF 
(Richardson's) method 
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Jc(ky=(I-U. D)X(k-l}HX C 



(15) 



or 
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*0fc>=a<(t-l>HX (c-A T A x(k-l)) 



x(k)=x(k-iyux A T [b-A x(lc-l)] 



16 



where k is the iteration index and a is the relaxation 
coefficient chosen such that all the eigenvalues of 
[I-aA^A] lie within the unit circle. Eq. 17 can be rewritten 
as 



(16) whether the system is overdetermined, underdeterrnined or 
undetermined. However, the solution may not be unique. 

It is evident from eq. 20 that, if the error corrections, 

(17) A r M4)]=0 (22) 

then, the solution converges after the iteration. 

When a system is overdetermined, or undetermined, i.e., 
the rank of A is less than M. The null space of A r is then not 
empty, and 
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x{z} = [/-[/-oo4 r A]r-q-%A*2>{z} 



(18) 



where z" 1 is the unit delay in terms of k and b{z} is a 15 
constant vector. Invoking the final value theorem, the final 
value of x is 



lim x(k) = 



lim[(l 
z=l 



(19) 



The result converges to the least squares solution of the 
normal equations as given in eq. 13. In tomographic image 
reconstruction, the matrix A is typically very large and 
sparse and consequently solving for x directly is unrealistic. 
The SIRT is a well-known method and uses a scheme 
identical to eq. 18 for solving an overdetermined system. 
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x(k)=rtk-iy¥ a 0 A T [b-A *(*-!)] 



c=M r 

x(k)=x(k-iyHX [c-$ A T A x(k-l)} 



The convergence properties of reconstruction algorithms 
expressed by eq. 20 and eq. 21 are the same. Both approach 
the least squares error in the projection space. However, 
when the system is overdetermined (M>N), the dimension of 
the back projected vector c is less than the dimension of the 
originally measured vector b. 

If numerical noise is not present, the SIRT expressed by 
eq. 20 guarantees convergence to a solution irrespective of 



there exists £?J_range of A, (23) 
where b - Ax(k) = e * 0, 
such that 
A T • e = 0, 

and A T [b - A*(£)] = 0. 

In this case, the rank of A is less than M and the residuals 
in the projection space may not converge to zero. However, 
the errors in the image space converge to the null space of 
the matrix A r . To emphasize the point further, the error 
corrections in the image space converge to zero even though 
the residuals in the projection space may not converge to 
zero. 

When M=N or M<N, the system is said to be exactly 
determined or underdeterrnined respectively. In these cases, 
there exists at least one k satisfying 



b-A *(*)=0 



(24) 
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where p is used to normalize the back projection defined in 
eq. 11 and eq. 12. When the system is overdetermined, the 
solution converges to the least squares solution. 

The SIRT can be considered as a time-invariant linear 
feedback system as shown in the block diagram given in 
FIG. 12. The matrix- vector multiplication associated with 
the matrix A performs the forward projection. The forward 
projection maps the reconstructed image into the projection 
space. This is equivalent to determining the measured pro- 40 
jections taken from the test object as described by eq. 9 and 
FIG. 10. The matrix-vector multiplication pA 7 ^ represents 
the back projection. The back projection maps the data from 
the projection space back to the image space. This operation 
is in the reverse direction of the forward projection and 45 
simply smears the data in the projection space back to the 
image space along each ray. The sparse matrix A is too large 
to be stored. The matrix-vector-multiplications associated 
with the back and forward projections are performed using 
eqs. 8 and 10 and weight factors, which are calculated on 50 
line (Wy's). Expressing SIRT using eq. 20, the residuals are 
obtained in the projection space as shown in the system 
block diagram in FIG. 12.3. From eq. 16, an alternative 
version of the original SIRT can be obtained by comparing 
the errors in the image space. 



When the rank of A is equal to M, there exist solutions for 
[b-A x] and the residuals in the projection space converge 
(20) 30 to zeros. Hence, the condition for convergence expressed in 
eq. 22 is still satisfied. However, when the system is 
underdeterrnined, i.e., M<N, the solution obtained using 
SIKT is not guaranteed to be unique since there are infinitely 
many values of x(k)'s satisfying eq. 24. An uncertain 
component in the null space of A consequently exists in the 
image space. The corresponding component in the solution 
vector, x, can be caused by computational errors accumu- 
lated during the iterative process. 
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Optoelectronic Implementation of SIRT 

In this section, an optoelectronic implementation of SIRT 
using SLM and CCD arrays is described in detail. The exact 
types of the SLM and CCD arrays are not specified in the 
proposed schemes. For example, the SLM arrays can be 
liquid crystal television displays and the CCD arrays can be 
solid state image detectors used in television cameras. 

As shown in FIG. 10, ray-sums of a projection are 
obtained from strip integrations of pixel gray levels along 
areas intercepted by grids of pixels on the image plane and 
paths of rays in the projection. A forward projection from the 
reconstructed image is a duplication of the operation per- 
formed on the original image except that it is taken from the 
reconstructed image rather than the object under test. In FIG. 
14, a 2-D SLM array is placed in front of a 1-D CCD array 
at an angle. The integrations are performed simultaneously 
by projecting the 2-D image pattern from the 2-D SLM array 
into the 1-D strip detector array. The 2-D SLM array 
represents the reconstructed cross-sectional image and the 
1-D CCD array collects the projected pattern of the image 
according to the geometrically intercepted fractional areas 
along each strip at the given angle. The operation on the 
SLM/CCD array pair merely duplicates the original projec- 
tion operation shown in FIG. 10 except that it is taken on the 
reconstructed image. The optoelectronic structure allows the 
matrix-vector-multiplication associated with the forward 
projection matrix at the q th angle, A q , to be performed by a 
single optical projection between the SLM/CCD pair. 



x„„<-x^,„+a A, 7 " [6,-A, (25) 2 " 
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A back projection, which is just an inverse of the forward of the optoelectronic devices have to be taken into account 

projection operation, smears the data from the 1-D projec- Several sources of degradation are discussed in the section, 

tion space back to the 2-D image space with the same source Q f degradation is caused by the finite 

intercepted i fractional areas between projection strips and ^ rf ^ deyices ^ ^ d ^ 

image pixels. In FIG. 13, a 1-D SLM array is placed in front 5 \ . \. * *. . ± A « * ^ , ~ 

of a 2-D CCD array, at an angle. The 1-D strip SLM array ««c s due to fact &at to « onl y ^ number of 

represents data in the projection space at the given angle and distinguishable gray levels which the optoelectronic device 

the 2-D CCD array represents the reconstructed image. After can generate. The second source of degradation is related to 

the light source is pulsed, the 1-D projected data is super- the nonlinear nature of the devices. The degradation in 

imposed on the previously reconstructed image in propor- lQ performance can also be caused by non-uniform distribution 

tion to the intercepted areas. The structure allows the matrix- 0 f light intensity incident on a plane from a point source. The 

vector-multiplication associated with the back projection _ " . . . , , 

matrix at theq rft angle, A q T , to be implemented using a single effec ? f random gyrated noise, which may be caused by 

optical projection between the SLM/CCD pair. the dark current of the CCD array, is also discussed. 

The SIKT algorithm of eq. 20 can also be expressed in The number of distinguishable gray levels associated with 

terms of projections taken at different angles: 15 CCD's is usually considered to be higher than those of 

repeat SLM's. Hence, only the effects from SLM arrays are dis- 

x OJd <—x ru;w cussed here. For a typical liquid crystal television display, 

for q=l to np do about two hundred distinguishable gray levels can be 

obtained. 

The effect of the finite dynamic range of the back pro- 
end of for jection SLM can be nnnimized by scaling the magnitudes of 
until done tne error tenns U P 311(1 i ^ ien seating the exposure period down 
where A T q , h q and K q are submatrices or subvector of A T , b correspondingly in steps 6 and 7 of the algorithm. In the 
and A corresponding to the q th angle respectively and np is 25 feedback configuration, the magnitudes of the error correc- 
the total number of projections from all the measured angles. tion terms decrease as the number of iterations increase. 
The structure is shown in FIG. 15. Forward projections are Therefore, the distortions caused by finite dynamic ranges of 
performed by the pair of arrays SLM2 and CCD2. Back the device can be minimized by scaling the magnitudes of 
projections are performed by the pair of arrays SLM1 and ^ ^ errors up to the full range of the back projection SLM 
CCDl. and then obtaining the proper error corrections by control- 
The iterative procedure using the structure shown in FIG. ^ the ^ interyal of me accordingly. The factors of 

15 can be described as follows: - . , . « £ « >4 . ^ 

* e «. scaling can be estimated in advance for each iteration. The 

procedure process for an iteration , . & . . , , . t . 

constant back, projection SLM and CCD operate at their maximum 

b • q=l to np, measured projection at the q th angle; 35 dynamic ranges, although magnitudes of error corrections 

variable keep decreasing as the number of iterations increase, 

x^: in CCDl, reconstructed image; Consequently, the number of distinguishable gray levels of 

x^: in SLM2, reconstructed image from the last iteration; a reconstructed image does not depend on the back projec- 

1. Postprocess the image obtained from the last iteration 40 tions any longer and instead depends on the number of 
(contents of CCDl) and then load it into SLM2 and CCDl distinguishable gray levels of the forward projection opera- 
as x old and the initial value of x new respectively. tions only. Therefore, the performance of the proposed 

2. for q=l to np do structures will rely on the 2-D SLM array used for forward 

3. Rotate SLM1 and CCD2 to the q* h projection angle; projections. The image in the forward projection SLM is 

4. Forward project x old , the reconstructed image from the 45 update(i on i y once during each iteration and consequently 

5. S errorTin projection space, [b — A x M ], by as fr ° m the perspective of overall system processing 
subtracting the reconstructed projection from the mea- s P eed * ^ less stongent response time requirement results 
sured projection. in more cnoices for selecting the device. 

6. Load the errors into SLM1 after scaling their magnitudes SLM's and CCD's are not perfectly linear devices. In 
up to the maximum dynamic range of the device. general, the degree of nonlinearity of SLM's is usually 

7. Scale the back projection exposure period down by the worse than those of CCD's. FIG. 16 shows the relationship 
same ratio as the scaling up factor used in step 6 and back between the input voltage and the transmitted intensity for 
project error corrections from the 1-D SLM1 into the 55 GaAs/AlGaAs -based CCD/MQW SLM's. In the block dia- 
reconstructed image in the 2-D CCDl. The result in gram of FIG. 17, block N(k)j5 and N(kV can be treated as 
CCDl is nonlinear functions of the back projection and forward 

projection data corresponding to SLM1 and SLM2 respec- 

jw=.*u + a i Aftbi-Aixou] (26) tively. The system can be described as 

" 60 

8 end of for «k>^ik-l)+A T N(k) B [b—A N(k) F Mk-i\=U-A T N{k) B A N(k) r ] 

9, end of the iteration. «M>** Wfc * (27) 

A scaling factor is used in steps 6 and 7 to eliminate the , _ a . , . _ 

effects of finite dynamic range. The factor will be described where and N ( k ^ m ^onal matrices reflecting the 

in the following section. 65 nonlinear functions. 

Since the proposed structure is essentially analog in If the system is stable, then, the final value of the 

nature, degradation in the performance caused by limitations reconstructed image is given by 
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lim x(k) - {ATNbANfY^A T N B b 

A'->oo 
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(28) 



= N^ 1 [A T N B .A]- 1 A T N B b 

where N B and Np correspond to the final value of N(k) B and 
N(k)„. 

In eq. 28, N 5 is neglected since its effect on the solution 
is secondary. Instead, we consider the effect of the nonlin- 
earity on the back projection first. The term [A r A]" 1 A r 
N £ b is equivalent to a least squares solution of a system 



N 2 Ax = N 2 

B B 



which minimizes (29) 



B 

The effect of nonlinearity on the back projection affects the 
distribution of the residuals. Therefore, the distortion asso- 
ciated with back projections can be neglected since the 
magnitude of the residuals are much smaller than the mag- 
nitude of the solution. 

The effect of the nonlinear nature of the forward projec- 
tion SLM cannot be neglected since it lies in the signal path. 
However, the effect is only to stretch the gray levels of 
reconstructed images. When a gray level of a pixel lies in the 
region of saturation of the SLM, the magnitude of the pixel 
is clamped at a constant value. Hence, the nonlinear char- 
acter of the SLM functions acts as a constraint on the values 
of the solution. The constraints on the values of the solution 
tends to make the system relatively more overdetermined. 
The net effect of this constraint is to ensure that the least 
squares solution is obtained in a reduced space. 

Another source of distortion is caused by the nonuniform 
distribution of the incident light intensity over SLM arrays. 
This can be caused by the Gaussian beam dispersion char- 
acter of lasers. Let T> B and Dp be NxN diagonal matrices 
representing the light intensity distributions of the back and 
forward projections respectively and let C(k) be an Nxl 
vector with identical elements representing the value of zero 
bias added during the back projection operation at the k* A 
iteration. The block diagram is shown in FIG. 17 and the 
system can be described as 
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xik)=x(k-l^-D B A T [b+C(k)-A D F x(k-l)]-A T C(k) 



(30) 
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where the zero bias C(k) is added to the data to be projected 
before back projecting and subtracted from the image after 
each iteration. The final image is given by 



K m 



x(k) = [DbATADfYIIDbA^ + IDb-IWCI 

= Dp\ATAY^ATb + Dp 1 [A T A\~ l rTp [Db - I\A T C 



(31) 



C- lim C(k). 



The term in eq. 31 associated with C is caused by zero 
bias and its effect can be reduced. As the magnitudes of zero 
10 bias is adjusted down as the reconstructed image converges 
to the final solution, the magnitude of C can also be reduced. 
If positive and negative valued error corrections are pro- 
jected separately, the zero bias term tends to cancel out The 
Dp -1 corrects the distortion caused by Dp at the forward 
projection and the distortion is canceled out if an output 
image is taken from the image used for forward projections. 
If the light intensity decreases as we move from the center 
to the edge of an image, the term associated with Dp -1 
causes the brightness of the reconstructed image to increase 
from the center to the edge of the reconstructed image 
directly obtained from the back projection CCD array. 

A dark voltage is generated in CCD cells as a result of 
build up caused by dark current. The effect of the dark 
voltage can be minimized by subtracting the average value 
from the CCD output after each iteration. Then, the residuals 
can be treated as random noise. 

If the system is overdetermined or exactly determined, 
then the effect of this random noise can be minimized by the 
feedback scheme of the SIRT algorithm. When a system is 
underdeteimined, i.e. the number of ray sums, M, is less than 
the number of pixels, N, the matrix [A r A] becomes singular. 
Components of the noise perpendicular to the null space of 
the matrix A will remain and accumulate in the image space. 
Consequently, convergence of the SIRT algorithm is not 
guaranteed. The problem of poor convergence is addressed 
hereinbelow where a new algorithm which guarantees con- 
vergence to a unique solution when the system is underde- 
termined is presented. 

In order to assure the validity of the scheme, the system 
was simulated. The effects of limited dynamic range, 
nonlinearity, and random additive noise associated with 
optoelectronic devices have been included in the simulation 
model. Results of the simulation are presented in this 
section. The Shepp and Logan phantom is used as a test 
cross-sectional image throughout the simulation exercise. In 
order to model the intercepted fractional areas of optical 
projections in the free space, the weight factors, w^-'s in eq. 
8, are calculated exactly in the simulations. The mapping 
errors, introduced as a result of the fact that the original 
measured projections are taken from continuous cross- 
sectional plane and the images are reconstructed on a 
discrete plane, are included. Projections are usually taken 
from an image with higher spatial resolution. The image is, 
60 however, reconstructed with lower spatial resolution. 

The Shepp and Logan phantom is a commonly used test 
cross-sectional image in tomographic image reconstruction. 
The phantom is specified in FIG. 18 and Table 3. Projections 
65 can be calculated directly from the given parameters of the 
ten ellipses or obtained by taking discrete Radon integra- 
tions from the image constructed using the given data. 
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TABLE 3 
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Parameters of the Shepp and Logan Phantom 



Center 


Major 


Minor 


Rotation 


Refractive 


Coordinate 


Axis 


Axis 


Angle 


Index 


(0, 0) 


0.92 


0.69 


90 


2.0 


(0, -0.0184) 


0.874 


0.6624 


90 


-0.98 


(0.22, 0) 


0.31 


0.11 


72 


-0.02 


(-0.22, 0) 


0.41 


0.16 


108 


-0.02 


(0, 0.35) 


0.25 


0.21 


90 


0.01 


(0, 0.1) 


0.046 


0.046 


0 


0.01 


(0,-0-1) 


0.046 


0.046 


0 


0.01 


(-0.08, -0.605) 


0.046 


0.023 


0 


0.01 


(0, -0.605) 


0.023 


0.023 


0 


0.01 


(0.06, -0.605) 


0.046 


0.023 


90 


0.01 



10 



15 



30 



35 



In simulations, an image is built first and then projections 
are taken from the image using eq. 8. The number of gray 
levels in the image is 200 as specified originally and a plot 
of a cross-section is shown in FIG. 19. However, the details 20 
of this phantom can not be observed if the gray levels are not 
stretched. The number of grey levels in the original phantom 
is 200, and the reconstructed image consists of real-valued 
gray levels in the range of 0-255. In order to enhance the 
contrast of the displayed image, gray levels in the range of 2 s 
126-134 were stretched after reconstruction to cover the full 
range of the display by multiplying by a factor of 32. FIG. 
20 shows a plot of the cross-section after contrast stretching. 

FIGS. 21a-d show reconstructed images as a function of 
the number of iterations. These images are reconstructed 
using the basic SIRT algorithm. Projections are obtained by 
strip integrating the 512x5 12 original cross-sectional image. 
The dimensions of reconstructed images are 128x128. The 
system is overdetermined (M>N) with the number of pro- 
jections set at 180. The width of ray-sums is the same as 
those of the pixels. Each projection covers the whole image 
plane. 

FIGS. 22a-d and FIGS. 23a-d show the results obtained 
by simulating SLM arrays with 1024 and 256 distinguish- 
able grey levels respectively. The finite number of gray 
levels is simulated by scaling the maximum pixel value to 
the maximum dynamic range and then quantizing properly. 
FIGS. 23a-d show that, when the number of distinguishable 
gray levels of the forward projection SLM array is reduced 
to 256, the reconstructed images show severe distortion. The 
results show that, for implementing the SIRT, the critical 
device parameter effecting the quality of reconstructed 
images is the number of distinguishable grey levels of the 
2-D forward projection SLM array. 

In this simulation, the nonlinear response curve in FIG. 16 
is approximated using a fifth order polynomial 



fa) = [y - 3.808249] x 202.3996 (32) 
where 

y = 0.00003205125 -0.003374 lz 4 + 0.063091c 3 -OAUOSz 2 + 
1.1307z + 2.46 
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2.0 10.0 and z = 5 £g +4.5 



The nonlinear function is employed during the back and 60 
forward projection operations separately. 

In this simulation, the system is underdetermined. The 
dimensions of the reconstructed image are 256x256. The 
number of projections are 120 and there are 121 ray-sums in 
each projection. Therefore, the ratio of M to N is 0.22 where 65 
M and N are dimensions of the matrix A in eq. 20. Uniformly 
distributed noise in the range of -0.5 to +0.5 is added to the 



reconstructed image after each iteration. The SNR is 48.27 
dB. FIG. 24 shows the plot of the residuals corresponding to 
the number of iterations. FIGS. 25a-d show the recon- 
structed image at different iterations. The result clearly 
begins to diverge after approximately one hundred itera- 
tions. 

The Simultaneous Iterative Image Reconstruction 
algorithm — SIRT has been discussed from the point of view 
of the method of least squares solutions for overdetermined 
systems. A scheme for implementing the algorithm is also 
presented. 

Optical computing offers the major advantages of massive 
parallelism and free space connectivity. These advantages 
can be exploited fully in tomographic image reconstruction. 
In the optical implementation of tomographic image recon- 
struction algorithm, the time consuming matrix-vector- 
multiplications associated with the back and forward pro- 
jection operations are replaced by parallel optical 
projections. The inaccuracy associated with interpolations 
and approximations linked to tomographic projection opera- 
tions are overcome since the intercepted fractional areas of 
free space optical projections are exact. In addition, iterative 
algorithms offer two major advantages over direct algo- 
rithms when implemented using optoelectronic devices. 
With iterative algorithms, the dynamic range of optoelec- 
tronic devices can be fully exploited. In addition, there is no 
need for filtering or convolution operations. The optoelec- 
tronic structures described in this section can be built using 
TV devices and mechanical image rotators that are com- 
mercially available. In the future, the performance of the 
proposed structures can be considerably enhanced by using 
higher speed SLM's and static image rotation devices. 

The presence of noise and optical distortion prevents 
optical implementations of direct reconstruction algorithms 
using optical Fourier Transforms with an acceptable level of 
performance. However, the accuracy of the reconstruction 
can be enhanced using the feedback scheme proposed in this 
section. In other words, a less accurate optical filter can be 
used with the optoelectronic implementation in order to 
speed up the reconstruction. 

In many applications, reconstruction of high spatial reso- 
lution images from fewer measured data is desired. The 
situation represents an underdetermined systems since the 
number of unknowns to be solved for, is more than the 
amount of measured data. It has been pointed out that, for an 
underdetermined system, the conventional SIRT algorithm 
does not guarantee unique convergence if numerical errors 
are present during computation. This problem is remedied in 
a new technique called Projection Iterative Reconstruction 
Technique which is described hereinbelow. An optoelec- 
tronic system for implementing the algorithm follows intro- 
duction to the novel approach. 

Projection Iterative Reconstruction Technique 

The Projection Iterative Reconstruction Technique (PIKT) 
is an iterative tomographic image reconstruction scheme for 
solving underdetermined systems. The PIRT guarantees 
convergence to a unique minimum-norm solution for an 
underdetermined system Since, the state matrix associated 
with the PIRT is symmetric and positive definite, it allows 
application of the conjugate gradient method for solving 
underdetermined systems without the need for imposing 
explicit constraints. 

In many applications of tomographic image 
reconstructions, images of higher spatial resolution are 
desired to be reconstructed from limited data. In all of these 
cases, the number of unknowns to be solved for, is more than 
the number of measurements. This results in an underdeter- 
mined system. 
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When a system is underdetermined, the SIRT algorithm The minimum-norm solution of an underdetermined sys- 

described hereinbefore converges to the minimum-norm tern can also be derived by minimizing (Vfc)x r x under the 

solution, if the initial value is selected properly and no equality constraint Ax=b. 
numerical error is introduced during computing. However, 

numerical errors are unavoidable in most situations. This is 5 d ( i ] (36) 

particularly true for the optoelectronic implementation pro- dx ( T ^-(Ax-Vfk j =0 

posed in the previous section. In this case, the result obtained T 

using SIRT will diverge as discussed hereinbefore and A*\^x = ° 

shown in FIG. 25. Consequently, there is interest in other AA T X=Ax = b 

algorithms where a unique solution is guaranteed when the 10 then 

system is underdetermined. ^ = [AA^b 

x = A 7 ^ 

The conjugate gradient method is an efficient acceleration 

technique for solving large and sparse linear systems and has = A^AA^b 
been used extensively in conjunction with the SIRT. Since, 

the conjugate gradient method is guaranteed to converge 15 where the Lagran ian multiplier X is equivalent to the state 

only when a system is symmetric and positive definite, it vector t in eq 33 
cannot be applied for solving underdetermined systems 

without imposing explicit constraints. Since the state matrix The minimum-norm solution of an underdetermined sys- 

associated with the PIKT is symmetric and positive definite, tern for the case of M=l, N=2 is shown in FIG. 26. 

the conjugate gradient method can be directly applied with- 20 When a system is underdetermined, i.e. M<N and all M 

out the need for constraints. rows G f A are linearly independent, the SIRT cannot guar- 

The PIKT is an iterative tomographic image reconstruc- antee convergence to a unique solution since the NxN matrix 

tion algorithm which can be considered as a counterpart of [A r A] is not positive definite anymore. However, in eq. 33, 

the SIRT. The SIRT leads to the least-squares solution of the tne MxM svrnme tric matrix [A A 7 ] is positive definite and 

normal equations for overdetermined systems by correcting . . - U1 . • j» * 

. i_ i .* . •* x- i r™. . therefore the state variable vector, i.e. the intermediate 

errors in the solution in the image space iteratively. This is , ^ < ~ . . , _ 

in contrast to the PIRT which yields the minimum-norm solution vector t, can be solved for lteratively. Consequently, 

solution of the normal equations for underdetermined sys- the minimum-norm solution of an imderdetermined system, 

terns by building up a state space vector in the projection Ax=b, can also be solved for by iteratively solving the 

space iteratively. 30 equation [A A^|t=b. The corresponding state and output 

Consider a system, A x=b, where M<N. If the M rows of equations can be written as 
the MxN matrix A are linearly independent, we have an 

underdetermined system and consequently there are an ^)#-i)4<xMA r <(M)] 
infinite number of solutions. However, the system has a 35 

unique minimum-norm solution. The solution of an under- and rtk)=A T t(k) (37) 
determined system which rmnimizes [fx|| 2 can be derived as 

follows where k is the iteration index and a is the relaxation 

coefficient, a has to be chosen such that the eigenvalues of 

If x—A?"t (33) 4Q [i-cc A A 7 ] lie within the unit circle. 

where t = [ttt 2 . . . t M ] T A block diagram summarizing the steps involved in PIRT 

then A[A T t] = b is shown in ^G. 27. The convergence of eq. 37 can be 

proved by recursively using eq. 34 and finding the limit as 

or t^lAA^b k-^~>. 
x = A T iAA T t 1 b 45 

«k)=t(k-l)+a[b-AA T \t(k- 1) (38) 

The matrix [A A^j is symmetric and positive definite and = ab + [/ - oAA 7 )^ - l) 

therefore a unique solution of t exists. = ctfr + {[/ - aAA^ab + [/ - vAA^ttjc - 2)} 

The fact that x represents the rmnimum-norm solution can = a{/ + [/ " aAA ' r ^ b + V ~ vAA 7 ) 2 ** - 2) 

be shown as follows. The solution of x is a linear combi- 50 f k-i l 1 

nation of the linearly independent rows of A, = « | ^ [/- aAA 7 }* j b + [I- aAA^O) J 



= a {/ - [/ - aAA 7 ]}- 1 {/+[/- aAA^b + [/ - oAA T \*t(0) 
= [AA T T X {/-[/- vAA r t\b + [/- aAA^tiO) 



x = A T t (34) 
= tiai T + tiai 7 + . . . + tMaM T 

where a, is the i th row of A. 55 When IP-ocAA 17 ]!^, and k->oo, then 

The solution vector x minimizes the Euclidian distance from k 

the origin because it is orthogonal to the null space of A. [i-aAA 1 ]*-*) (39) 



This can be shown as follows: 

Let ze sfr N and ze null space of A (35) 
then Az = 0 



60 



and 

(40) 



Urn t(k) — {AA T \~ l b 
a?z = 0, for all i,ie{ 1,2, . . . M} 

or * T = 0 65 for any finite t(0). 

consequently xlz where ± denotes orthogonality, i.e. xlnull The convergence can also be shown by taking 

space of A. z-transforms and invoking the final value theorem. 



5,654,820 



Kz) 



25 



Kz)2T l + a[b(z) - AA^z- 1 ] 



26 



[I-U-aAA^z- 1 ]- 1 - 



ab 



Km «k) = lim [(l-^Mz)] 

*-»oo Z— >1 

= [A4 T ) _1 fc 



Hence lim x(k) = A T lim 



A x=b, 

< 41 ) b Q represents the projections obtained from the initial image, 
b is the difference between the measured projections and the 
projections from the initial image, and x is the difference in 
5 the image space based on the information from the actual 
measurements and the initial image. The initial residual in 
the projection space is 



10 



= A T [AA T ]- 1 b 

This indicates that x(k) converges to the minimum-norm 
solution of the underdetermined system as shown in eq. 33. 
Since the state vector t is equivalent to the Lagrangian 
multiplier X in eq. 39, this algorithm can be considered as an 
iterative method for determining the Lagrangian multiplier 
and the optimum solution simultaneously. 

Alternatively, we can determine the minimum-norm solu- 
tion of the cross -sectional image of an underdetermined 
system using the PIRT algorithm by rewriting equation 37 as 



t(k)=t(k-iy+a [b-A x(k-l)] 
A T t(k) 



(42) 

where p is the coefficient used to normalize the back 
projections and t(k) is a state vector which is iteratively built 
up in the projection space. In this algorithm, the Lagrangian 
multiplier in the projection space and the minimum-norm 
solution in the image space are solved for simultaneously. 
The difference between the PIRT algorithm and the SIRT 
algorithm in respect of their computational procedure is that, 
with PIRT, there is an additional state vector in the projec- 
tion space and the error corrections are fed back to the 
projection space instead of in the image space. 
Consequently, the solution in the image space is uniquely 
determined by back projecting the reconstructed projections. 
In contrast the SIFT does not identify a state vector explic- 
itly and the error corrections are fed back in the image space. 
Consequently, in the case of PIRT, a unique minimum-norm 
solution is guaranteed when the system is underdetermined. 

In many applications, a priori information relating to the 
cross-sectional image can be obtained before applying itera- 
tive algorithms. An initial image can also be used in con- 
junction with the PIRT in order to improve the reconstruc- 
tion result as well as speed. When an initial image is given, 
the solution obtained using PIRT minimizes the Euclidian 
distance to the initial image in the image space. This is in 
contrast to SIRT which minimizes the residual in the pro- 
jection space. 

When the initial image is known, we only need to solve 
for an image that represents the difference corresponding to 
the actual measurements and the initial image. This differ- 
ence in the image space can be obtained from the initial 
residual in the projection space. The initial residual is found 
by projecting the initial image into the projection space and 
comparing it to actual measurements. Finally, the recon- 
structed image is obtained by superimposing the difference 
on the initial image. 

Let x Q denote the initial image. We rewrite A x=b as 



A [x 0 +x]=b 0 +B 



where 
A x 0 =b 0 



35 



40 



45 



50 



55 



(43) 



b = b-b 0 
— b —Axq 



(44) 



Therefore, only x, the difference in the image space, need to 
be solved for in the iterative reconstruction process. 



15 



and 



20 



m=f(k-ry^B-A z(k-i)] 

x(ky=A T f(k) 



(45) 



This algorithm minimizes ||x(.)|| 2 where 



25 



^.) = lim x(k\ 
k— >~ 



i.e., the Euclidean norm of the differences between the 
solution' and the initial image. 



30 



l|x-Xo|] 2 



(46) 



It may be noted that, it is not necessary for the initial image, 
Xq, to be the minimum-norm solution of A x 0 =b 0 . 

Constraints are very often used in conjunction with the 
SIRT and other iterative image reconstruction algorithms. 
Although a variety of constraints can be applied, only those 
applied to the values of solutions are discussed in this 
section. The constraints applied to the values of a solution 
confine the solution within a boundary that is usually based 
on an appropriate physical reasoning. Such constraints can 
include non-negative values or bounds for the pixel gray 
levels. Constraints imposed in conjunction with SIRT 
employed for solving an underdetermined system may not 
always lead to a unique solution. 

When constraints are applied in conjunction with the 
PIKT for solving an underdetermined system, the solution 
does not get affected if the unconstrained solution is within 
the solution boundary. If the unconstrained minimum-norm 
solution is outside of the boundary, the solution converges to 
a minimum-norm solution in a reduced solution space 
provided that the minima exists within the reduced space 
and the constrained space is convex. In this case, the solution 
x can be expressed as [x v T !x c Y, where x v is an N^l vector, 
x,, is an N c xl vector and N V +N C =N. x c is a vector of constant 
values which corresponds to the part of the solution which 
is constrained by the boundary. x v denotes the rest of 
variables to be solved. The problem is equivalent to solving 
a system with reduced dimensions 
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[A v \A e ] 



A v x=b-A c x=c 



(47) 



(47) 



65 where [A V !A C ]=A and c=h-A c x c . 

The dimension of the reduced system, A v x v =c, is MxN v 
where N V <N. The solution converges to the unique 



27 



minimum-norm solution in the reduced space provided that 
the system is still underdetermined. 

When more and more pixels are constrained to take 
constant values, the dimensions of the solution space may 
become smaller than the dimension of the measurements in 
the projection space. In this case, the originally underdeter- 
mined system is turned into an overdetermined system 
subject to the constraints. It is also possible that, after some 
pixels are constrained to take constant values, the number of 
linearly independent columns in the matrix A is less than M 
despite the fact that more than M pixels need to be solved 
for. In the latter case, the system become an undetermined 
system. In both of the above cases, the matrix [A v A y T ] is not 
positive definite. 

The state vector t in the projection space of the PIRT 
expressed by eq. 42 is only guaranteed to converge to a 
unique solution when the system is underdetermined or 
when the matrix A is square and full rank. The residual in the 
projection space will converge to zeros only when the rank 
of A is M. When the rank of A is less than M, the state 
variable, t, in the projection space will diverge as shown 
below. 

Consider an overdetermined system, where the measured 
vector b can be decomposed into two orthogonal compo- 
nents 
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and 



t ± (ky=t ± (k-l}+[b ± -A A T r x (*-l)] 



where 

t,|(k) e range of A 
tj^k) e null space of A T 

for r fl (0)^ x (0)=0 



Consider a measurement b± in the null space of A r , b±?K) 
and A T b±=0. The vector t_L therefore linearly increases in 
proportion to the iteration number k: 
since A T t x (k-l)=0 



= t ± (k-i) + [b ± -AA T t ± (k-iy\ 
= t±(k-i)+b± 

= t ± (k-2) + 2b ± 



forfj.(0) = 0 



lim \\t ± m 2 = 



This indicates that the state vector t has a tendency to diverge 
linearly as a function of the iteration number. 



10 



(48) 
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let b±e. sft M with 

b± null space of A 7 " 

and 

6||€ sft Af with 
b\\ e range of A 

then 

b = aix\ + azX2 + . . . 4- afpCN = b± 

= b\\ + b± 
for b\\ _L b± 



where M>N, and a y is the j th linearly independent column of 40 
A. Under these Conditions, the linear system can be decom- 
posed into two independent systems 
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(49) 
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(50) 



(51) 



Hm t(k)= lim [t\\(k) + t ± (k)] = <~. 



15 



20 



When the system is overdetermined, the solution obtained 
using PIRT in the image space converges to the least squares 
solution despite divergence in the projection space. 

Let x(.) denote the final converged solution and x(.) y , j=l, 
2, . . . , N, be the } th element of x(.), then 

b r^i x x +a 2 x 2 +. . . -¥a N x N 



since b,-A x(.)=0 
hence 

25 a x [^-^(.)i]+^2 [^-*0 2 ] + ■ - ■ +*n [**-*O*]=0 

and all a^'s are linearly independent 
therefore, x y — x(.)/=0, for all j=l to N. 
Then the residual is 

b ± =b-A x(.) 



(52) 



(53) 



The solution is equivalent to the one obtained by minimizing 
l|b-Ax(.)|| 2 . 

It has been shown that for the convergence of PIRT in the 
projection space, the system has to be underdetermined or 
determined. In this section, an important restriction on the 
projection geometry of the PIRT is introduced. It should be 
kept in mind that whenever the PIRT is used, none of the 
projections should cover the entire support region of the 
solution in the image space. 

In the case where M^N for an MxN matrix A, , a set of 
projections should be chosen so that all columns of the 
matrix A are linearly independent in order to avoid having an 
indeterminate system. It can be shown that, if any two 
projections cover the entire support region of the solution x 
on the image plane as shown in FIG. 28, then the columns 
of A will become linearly dependent. The system then 
becomes ^determinate. 

For an MxN matrix A, M^N, if a row of A can be 
expressed as a linear combination of other rows of A, i.e. for 
any a, e<R 



Via? 



(54) 



N} 



where &f is the i th row of A, i, p e {1,2, 
then the system is undeterminate. 

If we assume that there are projections at angle 8^ and B g 
which cover the entire support region of the variable vector 
x on the image plane, then both projections have the same 
total mass, Le. 



i P eS p 



-. X 



N 



(55) 



65 where and denote sets of indices of ray-sums corre- 
sponding to the projections at 0 p and 0^ respectively. 
Then 
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Therefore, the system is indeterminate since & p T is a linear 
combination of the other rows of A. 

If the solution x is constrained to be strictly in the convex 
support region on the image plane, defined by 



dk+i T A <**=0> for all k_ 



(66) 



Eq. 63 is derived by substituting eqs. 64, 63, 62, 60 and 59 
into eq. 66. (^represents the optimum step size in the linear 
search. The update equation, eq. 60, is derived from 



£ aj 



Z i 
i p eS p 



(56) 



dF(x k+ j) 



= 0 



(67) 



10 



i€S p vjS q 



a? 



the support regions of projections as shown in FIG. 29, then 
the system becomes undetermined. This is a consequence of 
the fact that every projection covers the entire support region 
of the solution x. In this case, constraints on the support 
region of solutions need to be somehow relaxed in order to 
prevent divergence of the state vector, t. 

Conjugate Gradient Acceleration of SIRT & PIRT 

The conjugate gradient method has been extensively used 
as an acceleration technique for solving linear systems. 
Unlike other iterative methods which may involve an infinite 
number of iterations, the conjugate gradient method, com- 
putes the exact solution in a finite number of iterations. 
However, the convergence is guaranteed only when the 
system is symmetric and positive definite. 

Consider the problem where we wish to minimize a 
quadratic function 



15 



20 
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where the function, F(x), is defined in eq. 57. 

Since, all the search directions are conjugate to all pre- 
vious search directions, the irdnimization along a given 
direction implies a minimization along all previous search 
directions also. Hence, for an N dimension problem, only N 
steps are needed to arrive at the solution provided there is no 
numerical error. If the system has M distinct eigenvalues, the 
conjugate gradient method converges in M steps rather than 
N steps. Examples of a two dimensional case are given in 
FIGS. 30 and 31. In the example, the conjugate gradient 
method arrives at the solution in at most two steps whereas 
the steepest descend method requires an infinite number of 
steps. 

In order to guarantee the convergence of the conjugate 
gradient method, the matrix A has to be symmetric and 
positive definite. Otherwise, the relations given by eqs. 58, 
65, 66 and 67 may not be valid. 

When a system is overdetermined, the matrix [A A 7 ] in 
eq. 13 is symmetric and positive definite. Hence, the con- 
vergence of the conjugate gradient method is guaranteed. In 
this case, the descent direction p* in eq. 59 is updated by 

p*-A r Ad* (68) 



F{xy=Y2x T A x-b T x 

The gradient of F(x) is 



(57) 



Then the iterative procedure is 
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Pk 



(58) 



. A T Ad k 
Pk T dk 



dk T Pk 



= Ax-b 

where, r(x) is also the residual of A x— b. When r(x) is zero, 
then the corresponding x is the solution of the linear system 
A x=b. 

The conjugate gradient method of solving the equation A 
x=b involves the use of following iterative procedure 



40 



r k+l * 



45 



d k T n 



Pk T d k 



d * T Pk 



X M «= x k + V.*d k 
r k+i «= r k - Cf^p k 



Pk T dk n T r k 
dk+j <= r k+i + PiA 

^59) The corresponding initial values can be chosen as Xo=0 and 
dQ=r 0 =A T b. Eq. 69 is the commonly used SIRT-CG type 
(60) algorithm in image reconstruction. 

When a system is underdetermined, the matrix [A A 7 ] in 
50 eq. 33 is symmetric and positive definite. The convergence 
of the conjugate gradient method is, once again, guaranteed 
In this case, eq. 59 is replaced by 



(61) 
(62) 

(63) 



Pk T dk 

4t+7 <= r k+j + Pk<4 



p*<-AA r cU 



(70) 



(64) 



The initial values can be chosen as x o =0, do=b and r 0 =b. In 
the iterative procedure, p* and a* are the descent direction 
and descent step respectively, d^. and represent the search 
direction and the search step. 

In the iterative procedure described by eqs. 59 to 64, the 
search direction d k is conjugate to all the previous search 
directions, i.e. 
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The state vector can be solved using the following procedure 

Pk *= AA T d k (71) 

dk T n n T r k 

or- 



Pk T d k 



dk T Pk 



dk+i T A dj=f), i=0, ...,k. 



(65) 



65 



*k+i •*= *k + 

r k + J OfcP* 

Q -pk T rk+i rf+in+i 
P* = ~ or- 



This is accomplished by choosing the search direction step 
size such that 



Pk T dk 



n T r k 
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31 

-continued 



where the corresponding initial values can be chosen as 1^=0, 
do=b and r 0 =b and the final solution is x=A T t. Eq. 71 
represents the PIKT-CG algorithm for tomographic image 
reconstruction. It can also be solved without building up the 
state vector 



Pk 



*= AA T d k 



(72) 



dk T n n T r k 
— — — or- 



d k T Pk 



r k+l • 



Pk T d k 



dk+i «- r k+j + Pk<4 

where the corresponding initial values can be chosen as 
x o =0, d Q -b and r c =b. The iterative procedure expressed in 
eq. 72 is similar to the algorithm used for solving unsym- 
metric systems. 

For large and sparse systems, it is difficult to obtain the 
matrix-inatrix product [A T A] or [A A 7 ]. In this case, the 
conjugate gradient method can be used directly on the 
matrices A and A r sequentially using matrix-vector- 
multiplications as shown in eq. 69. 

In the case of tomographic image reconstruction, the 
SIRT-CG algorithm was developed as an alternative to SIRT 
in order to improve convergence for both overdetermined 
and underdetermined systems. 

When a system is overdetermined, the matrix [A r A] is 
symmetric and positive definite. Hence, the conjugate gra- 
dient method can be applied in conjunction with the SIRT 
type algorithm with guaranteed convergence. Unfortunately, 
when the system is underdetermined, the matrix [A T A] 
associated with the SIRT is not positive definite anymore. 
Hence, the acceleration technique does not lead to conver- 
gence if any errors are accumulated during the iterative 
process. However, with the PIRT type algorithm, the state 
space matrix [A A 7 ] is symmetric and positive definite when 
the system is underdetermined. Consequently, the PIRT-CG 
methods can be applied for solving underdetermined sys- 
tems. Convergence of such algorithms is then guaranteed 45 
even if explicit constraints are not applied. 

Since the residual in the algorithm expressed in eqs. 69, 
71 or 72 is calculated recursively instead of obtaining from 
the current state t(k) (or x(k)) and the original input b, 
decoupling between the system and the original input can be 
caused by numerical errors or noise associated with opto- 
electronic devices. In order to avoid the decoupling, the 
residual x k+1 can be updated using the true residual, [b-A 
x k+1 ], instead of using the recursive approach. The PIRT-CG 
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algorithm can be rewritten as follows 
P*<-[A A 7 ) d k 
v*k=r k T rflk' Pk 

x* + i«-A T t Jt1 _ 1 
Jjfc+i<-b-A 
$k=r k+1 T r M /r k T r k 
dk+i*-*k + i+$k d * 



where t o =0, do=b and r G =b. 

The nratrix-vector-multipHcations associated with the 
PIRT-CG can be performed using optical projections. This 
allows implementation of the algorithm using a high speed 
optoelectronic structure. Implementation details are pre- 
sented in the next section. 

The algorithm of eq. 42 can also be expressed in terms of 
projections from different angles: 

repeat 



15 



20 



25 



30 



35 



40 



for q=l to np do 

t^+ct [ V A 9 x„ w )] (74) 

end of for 
until done 

where t^, A q T y b q and A q are submatrices or subvectors of t, 
A r , b and A respectively corresponding to the q th angle. The 
total number of projections from all measured angles is 
represented by np. 

Implementation of PIRT 

The proposed optoelectronic system for implementing a 
PIRT underdetermined system is shown in FIG. 32. The 
structure is similar to that for the SIRT shown in FIG. 14. 
The implementation of the forward and back projection 
operations were explained hereinbefore. The complexity of 
operations is about the same as SIRT but extra storage is 
needed for the state vector t in the projection space. In 
addition, the iterations are not carried out in the image space. 
Instead, a new image is reconstructed using only the updated 
state vector in the projection space. 

The iterative procedure which can be implemented using 
the structure shown in FIG. 32 can be described as follows: 
procedure 

constant 

b q : q=l to np, measured projection at the q th angle; 
variable 

t q \ q=l to np, state vector at the q th angle; 
x^: in CCD1, reconstructed image; 
x ald : in SLM2, reconstructed image from the previous 
iteration; 

1. Update x old in SLM2 with x Tiew in CCD1; reset x new in 
CCD1 to 0; 

2. for q=l to np do 



(73) 



so 3. Rotate to % 

4. Forward project x old , the reconstructed image from the 
2-D array SLM2 into the 1-D array CCD2, thereby 
generating [A q x old \. 

5. Update t q by {t^+tb^-A^ x old \}\ load t q into the 1-D array 
55 SLM1; 

6. Back project the state vector in the 1-D SLM1 into the 2-D 
CCD1. The result in CCD1 is a superimposed set of back 
projections: 



60 



i=l 



Ai T ti 



7. end of for 

8. end of the iteration. 

65 The conjugate gradient method can also be implemented 
using the structure shown in FIG. 32. In the accelerated 
algorithm described by eq. 73, there are two matrix-vector- 
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multiplications associated with the matrix [A A 7 ]. Using the 
optoelectronic structure, these matrix-vector-multiplications 
are performed by passing the vectors d and t through the 
back projection and the forward projection operators sepa- 
rately. The reconstructed image is obtained without addi- 5 
tional cost since the output vector x represents the interme- 
diate result obtained after back projecting t The rest of the 
calculations in the algorithm are simple vector-vector or 
scalar operations and therefore, do not require any special 
structures. 10 

The simulation results presented in this section are in 
agreement with the analytical properties with respect to 
PIRT. The performance of the PIRT has been verified in 
respect of convergence properties both for the basic as well 
as the modified version employing the conjugate gradient 15 
algorithm. 

FIGS. 33a-d show the reconstructed image obtained 
using the basic PIRT while FIG. 34 presents plots of 
Euclidian norms of the image and the residual as a function 
of the number of iterations. The residual is the distance 20 
between the reconstructed image and the original image. The 
reconstructed images are 128x128 and consequently N=16, 
384. The number of projections are 120 and there are 121 
ray-sums in each projection. Therefore, M= 14,520. The ratio 
of M to N for the underdetermined system is 0.886. In order 25 
to study the convergence characteristics, constraints have 
not been imposed. 

FIG. 35 and FIG. 34 show results obtained using the 
conjugate gradient method. FIGS. 35a-d show the recon- 
structed image while FIG. 34 gives the plots of the Euclidian 30 
norms of the solution and the residual as a function of the 
number of iterations. The corresponding values of the origi- 
nal images, the reconstructed images and the residuals 
obtained using the two methods are listed in Table 4. The 
images reconstructed using the basic PIRT are relatively 35 
smooth compared to those obtained using the conjugate 
gradient method. However, the residual of the conjugate 
gradient method at the $ th iteration (16.06363) is approxi- 
mately the same as that obtained using the basic method at 
the S2 th iteration ( 16.07527). The residual obtained using the 40 
basic PERT at the 8** iteration is 34.65206. 



guarantees that the solution converges to the unique 
minimum-norm solution in the image space. With a priori 
information, this algorithm minimizes the Euclidian distance 
between the reconstructed image and the initial image in the 
image space. When constraints are applied and the 
minimum-norm solution is outside the constraint boundary, 
the solution converges to a unique minimum-norm solution 
in the reduced solution space. It also allows the application 
of acceleration techniques, such as conjugate gradient 
method, directly to the underdetermined system without the 
need for imposing explicit constraints. 

When the system is underdetermined, the PIRT guaran- 
tees convergence in contrast to the SIRT which fails to 
converge. 

Iterative Filtered Back Projection — IFBP 

In this section, the Iterative Filtered Back Projection 
(IFBP) technique is presented. This algorithm is based on 
the strategy of iteratively applying the Filtered Back Pro- 
jection (FBP) method to approach the weighted minimum 
least square error in the frequency domain. In addition, a 
scheme for implementing the system using optoelectronic 
devices is given. 

The Filtered Back Projection (FBP) is considered as a 
direct method for tomographic image reconstruction. The 
FBP and the corresponding convolution algorithms are com- 
monly used in medical and industrial applications because of 
their lower computational demands. In some situations, 
reconstructed images can be distorted due to several factors, 
such as those introduced by aliasing from undersampling, 
finite filter bandwidth, limited views of projections, finite 
dynamic range of optoelectronic devices, and especially, 
speckle noise associated with optical transforms. The itera- 
tive approach introduced in this section can be used to 
improve the quality of reconstructed images by applying 
error corrections repeatedly. It will be shown that the errors 
introduced by the FBP procedures can be minimized in 
terms of the weighted mean square error in the frequency 
domain. Although FBP procedures can be considered as a 
method for estimating the inverse of the Radon transform, 
the accuracy is poor when the inverse is obtained using 



TABLE 4 



Residuals of Reconstructed Images Using PIRT and PIRT-CG. 
iteration 



algorithm 0 1 2 4 8 16 32 64 82 



PIRT 102.54 61.16 49.68 41.57 34.62 28.40 22.53 17.54 16.08 
PIRT-CG 102.54 57.32 38.97 26.51 16.06 11.41 10.35 9.80 9.58 



In this simulation, the system is severely underdeter- 
mined. The dimensions of the reconstructed image are 55 
256x256. The number of projections are 120 and there are 
121 ray-sums in each projection. Therefore, the ratio of the 
M to N is 0.222 where M and N are dimensions of the matrix 
A in eq. 42. Uniformly distributed noise in the range of —0.5 
to +0.5 is added to the reconstructed image after each 60 
iteration. The resulting SNR is 48.27 dB. FIGS. 36a-d show 
clearly that the PDRT converges monotonously even with the 
additive noise. In contrast, the SIRT algorithm fails to 
converge as shown in FIG. 25. 

The tomographic image reconstruction algorithm — PIRT 65 
is obviously superior to SIRT in respect of its convergence 
properties for underdetermined systems. This algorithm 



optical transforms. However, optical transforms can be used 
to accelerate the iterative reconstruction procedure based on 
the method proposed hereinbefore. 

The IFBP method is presented after a brief review of the 
FBP. The properties of the optical Fourier Transform are 
described first. An optoelectronic complementary filter 
structure is then introduced, where the bipolar radius filter of 
the FBP is implemented using a low pass optical comple- 
mentary filter. The optoelectronic implementation of the 
IFBP is presented and simulation results are given. 

The Fourier Slice Theorem can be stated as follows: The 
Fourier transform of a parallel projection of an image f(x,y) 
taken at angle 0 is equal to the 1-D slice of the two- 
dimensional transform, F(u, v), subtending an angle 0 with 



35 



the u-axis. In other words, the Fourier transform of P(6,p) is 
identically equal to F(u, v) along the axis u* in FIG. 37. 

In FIG. 37<s, a projection is taken on the x-y plane of a 
cross- sectional image by line integrating f(x, y) in a direc- 
tion perpendicular to the x' axis where the angle between the 
x' axis and the x axis is 8. As shown in FIG. 37&, the Fourier 
Transform of this projection can be obtained by computing 
F(x, y) along u' where the angle between axis u' and axis u 
is also 0. Consequently, the entire frequency plane can be 
established by taking projections at various values of 0 from 
0 to 71 and then transferring the information into the fre- 
quency domain. 

A cross-sectional image can be reconstructed using the 
2-D Inverse Fourier Transform to map the information from 
the u-v plane to the x-y plane as shown in eq. 75. When the 
Inverse Fourier Transform is expressed in a polar raster 
coordinate system, an additional term, p representing the 
radius in the frequency plane, is introduced in the integral. 
Since the integrations are separable, the first integration 
within the brackets is equivalent to the filtering at a given 
angle 0 using a filter whose kernel is [pi. 
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where, the matrix A T is the back projection operator matrix, 
the vector b contains the measured projections, the matrices 
F and F I2VV are linear operators representing Fourier and 
Inverse Fourier Transform operations respectively, the 
diagonal matrix R is the kernel of the radius filter and np is 
the total number of projections. The matrices and vectors 
with the subscript q correspond to the q* h projection angle. 

When the system is overdetermined, the FBP procedure in 
eq. 77 can be used iteratively to minimize the square errors 
in the frequency domain. An iterative procedure based on the 
basic Richardson (RF) method is given by the following 
difference equation 



x(ky=x(k-l)+a A T F lNV R F [b-A x(k-l)] 



(78) 



where k is the index of iterations, a is the relaxation factor, 
matrices A and A T are the forward projection operator matrix 
and back projection operator matrix as defined in eqs. 9 and 
11. FIG. 39 shows the block diagram for the iterative filtered 
back projection algorithm. 

The IFBP expressed in eq. 78 can be considered as a 
time-invariant linear multi-variable feedback system with 
step input. The reconstructed image is the steady-state 
output which can be found using the z-transform and invok- 
ing the final value theorem. Using z-transforms, the differ- 
ence equation in eq. 78 can be written as 



Q(Q,r)dQ 



where r=x cos 0-fy sin 0, (0, p) and (0, r) are the radii and 35 
the angles along the polar raster scans corresponding to the 
u-v and x-y planes respectively. In eq. 75, f(x, y) is the value 
of the image at the coordinate (x, y). F(u, v) is the Fourier 
representation of f(x, y) on the u— v plane. P(0, p) is the 
Fourier presentation of a projection at angle 0 and radius p, 40 
and is equal to F(u, v) on the 0-p plane. Q(0, r) is the filtered 
projection on the 0-r plane. As shown in eq. 75, the filtered 
back projection Q(0, r) at a given 0 can be obtained by 
computing only the 1-D Fourier Transform and filtering it 
using a 1-D filter kernel shown in FIG. 38. 45 

The representation of eq. 75 in discrete form is the 
Filtered Back Projection (FBP) method which can be 
expressed as 



*M = xtoz- 1 + aA T F INV RF[b\z\ - Ax{z\z~^ 
= [/-[/- uA T F ]NV RFA Jz~ [ ]aA T F fNV RFb - 
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Obviously, when the matrix [A r F /w RF A] is non-singular 
and a is chosen such that all eigenvalues of [I-a A T F INV R 
F A] are within the unit circle, eq. 79 converges to 



(80) 



lim 



lim x[z] - 
z->l 



= [A T F !NV RFAY l A T Fj NV RFb 

Since the Inverse Fourier Transform matrix F INV is the 
conjugate transpose of the Fourier Transform matrix F, i.e., 
F /Arv ^=F*, and the diagonal matrix R can be decomposed as 
R=(R) 7 R, the algorithm minimizes 



Ti Po (76) 50 

X(iuh) = 2 X 7>(e,p)lpteP*Pr 

e=o p=-po 

= 2 0(8,r) 

where 1=^ cos 9+i 2 sin 0, i t and i 2 are coordinates in the 
discrete image plane, and x(i l9 i 2 ) represents the grey level 
of the (i l9 i^) th pixel. In eq. 76, the 1-D filtering operation can 
be performed separately and then overlapped into the image 
plane. Since all the indices in eq. 76 are in the discrete form, 
interpolations are usually used. 

It has been shown in FIG. 10 that projections can also be 
represented by means of weighting factor w(i, j)'s which are 
the fractional areas intercepted by ray- sums of projections 
and pixels of the image. Therefore, the discrete form of the 
filtered back projection method of eq. 76 can also be 
expressed in terms of the back projection matrix, A T . Using 
linear operator matrices, the FBP procedure can be 



1 1 1 

||fl 2 FAx- R 2 Fbfo= \\{R 2 F){Ax- 



(81) 



where ||.|| 2 denotes the Euclidean norm of a vector. Eq. 81 
shows that the solution of IFBP is obtained by minimizing 
the least square error of the filtered projections in the 
55 frequency domain with a square root radius filter. 

The diagonal matrix R in eq. 78 represent the radius filter 
kernel as shown in FIG. 38. Unfortunately, the zero and near 
zero values in the radius filter operator inatrix, reduce the 
rank of the system. For an MxN system A x=b, where M and 
N are the column and row dimensions of the matrix A 
respectively, the dimensions of R are MxM. When the 
number of nonzero entries in the diagonal matrix R is less 
than N, the system is singular and convergence is not 
guaranteed. In order to have enough nonzero entries in R, the 
system has to be overdetermined, i.e., the number of pro- 
jection data should be larger than the number of pixels in the 
reconstructed image. 



5,654,820 



37 



As described hereinafter, the conjugate gradient method 
can be used as an acceleration technique for solving over- 
determined and underdetermined systems. When the system 
is not singular, the matrix [A r F r (R) r R FA] in the IFBP 
method is positive and symmetric definite. Therefore, the 5 
convergence of the IFBP can be accelerated using the 
conjugate gradient method A unique solution is guaranteed. 

Optical Fourier Transform Using a Two Lens 

System io 

In this section, the optical Fourier Transform using a two 
lens system is briefly described. An optical complementary 
filtering structure implementing the radius filter of the FBP 
is then introduced. 

Consider a two lens coherent optical system as shown in 
FIG. 40, where F x and F 2 are the focal lengths of the lenses 
respectively. If an image, f(x, y), is placed at the front focal 
plane of the first lens, its Fourier Transform in the spatial 
frequency domain, g(u,v), is obtained at the back focal 
plane: 



C\ re -frit 



""I dxdy 

- ° l F [ U V 1 
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where c t is a constant and 1 is the wave length of the 
monochromatic light. The image at the back focal plane of 
the second lens is a "flipped" Inverse Fourier Transform of 30 
g(u,v): 
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where D is the diameter of the aperture. 

For incoherent optical systems, the spatial frequency 
spectrum of the image can be obtained by applying the 
complex weighting factor H(p) to the spatial frequency 
spectrum of the object intensity. The function H(p) is known 
as the optical transfer function (OTF) and the modulus !HI is 
the modulation transfer function (MTF). 

The OTF for the aberration-free lens system of a circular 
aperture can be found to be 
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ff(p) = 
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(85) 
P = 2po 
otherwise 



where p Q is given in eq. 80. The spectrum of H(p) is shown 
in FIG. 42. 

The 1-D radius filter in the filtered back projection 
method has a frequency response given by 
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A(p) = 



P 
Po 



0 ^ p ^ p 0 



(86) 



where p 0 is the cut-off frequency. The frequency response is 
shown in FIG. 43. The filtered output is 



g(ry=IFT{FT{f(r)}.h(p)} 



(87) 



where g(r) and f(r) are the output and the input of the filter, 
respectively. 



A(V,y) = Jfe(«, v)e^ K XF 2 dudv 



(83) 
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where c 2 and c 3 are constants and F^/Ft is the magnification 
coefficient. 

If a spatial frequency mask is placed between the back 
focal plane of the first lens and the front focal plane of the 
second lens, h(x\ y*) represents a filtered version of the 
original image. 

In coherent optical systems, the representations of f(x, y), 
g(u, v), h(x f , y') and the filter mask are in terms of the 
magnitude of the light. In contrast, in an incoherent optical 
system, the intensity is used instead of the magnitude. In 
optical information processing, coherent optical Fourier 
filtering usually suffers from speckle noise caused by the 
granularity of imaging devices and dust. This is considered 
to be one of the major obstacles in optoelectronic imple- 
mentations of tomography. 

The finite passband in the Coherent Transfer Function 
(CTF), within which the coherent imaging system exists, 
passes all frequency components without amplitude and 
phase distortion. At the boundary of this pass band, the 
frequency response drops to zero, implying that the fre- 
quency components outside the passband are completely 
attenuated. The frequency response of a circular aperture is 
shown in FIG. 41. The cutoff frequency in terms of the 
image spatial frequency is 



Since h(0)=0, the zero frequency component should be 
eliminated completely and the dynamic range of the outputs 

45 must be bipolar. However, in an optoelectronic system, 
photo detectors can only detect intensity of light, and 
consequently both input signal and output signal should be 
biased to non-negative values. Therefore, the optical filter is 
required to be able to separate the zero frequency component 

50 of the bias and the zero frequency component of the signal. 
The zero frequency component of the bias should be passed 
while suppressing the zero frequency component of the 
signal. These requirements are contradictory. In order to 
circumvent the conflicting requirements, an optical comple- 

55 mentary filter is proposed. The complementary radius filter 
decomposes the original radius filter into two components, 
one is the filtered complementary function of the input while 
the other is the original input. The output of radius filter is 
obtained by superimposing the filtered signal on the original 

60 input The frequency response of the low pass complemen- 
tary filtering channel is given by 

Hp) = rf( P )-A(p) 
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where d(p) = 
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-continued 
Ipl ^ po 
otherwise. 
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The frequency response is shown in FIGS. 43a-b. The 5 
output of the filter can be described as follows 

gXr) = [ "T /FT { FT | -§■ } 'HP) j (89> 

= 4j- + IFT{FT{f(r)}} + IFT | FT j J • h(p) J 

/Fr{Fr^r)} ■ [-A(p)]> 
= c + IFT{FT{f(r)}'[<i(p)-Kpm 
= c + /Fr{FT{/tr)}-fc(p)} 
= c + g(r) 

where g'(r) is the output of the complementary filter and c/2 
is the value of the zero frequency component bias light 
intensity applied to both the filtered and unfiltered channels. 20 
The terms [c/2-ff(r)] and [c/2-f(r)] in eq. 89 represent biased 
nonnegative inputs applied to the optical filter. This filter can 
be implemented by a two lens system with symmetrical 
positive and negative response SLM's. FIGS. 44a-c show 
I/O response curves of a pair of symmetric positive and 25 
negative SLM's. 

The dynamic range of the filtered output is c±c/2 repre- 
senting about 33% loss in the dynamic range. However, the 
loss can be compensated for in the form of a zero frequency 
component bias in the input to the next optically addressed 30 
SLM. 

Using the complementary radius filter structure, the high 
pass spatial frequency response of the radius filter is 
replaced by a low pass response characteristic. In addition, 
the proposed feedback scheme offers increased tolerance to 35 
distortion introduced at the filter stage. The complementary 
radius can be implemented using incoherent optical systems, 
such as, the Ronchi pupil and optical transfer function (OTF) 
synthesis techniques. 

The iterative procedure of the IFBP in eq. 78 can be 40 
rewritten as a summation in terms of operations of each 
projection angle. Let q be an index of projections and q=l, 
2, . . . , np, and np be the number of angles, then 
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where u(k) is inserted as a scaling factor used to keep 
optoelectronic devices working at their full dynamic range 50 
and a is a relaxation coefficient. Variables in eq. 90 with 
subscript q represent the corresponding submatrices and 
subvector at the q th projection angle respectively. 

FIG. 45 shows an optoelectronic structure for implement- 
ing the IFBP algorithm. There are three major functional 55 
blocks in the structure which implementing the back and 
projection, the forward projection, and the optical radius 
filter. Implementations of forward and back projections are 
similar to those used for implementing SIRT and described 
hereinbefore. Amplifiers are inserted prior to optical filter- 60 
ing. Magnitudes of the error terms are scaled by the factor 
u(k) up to the full device range before filtering is performed. 
During the back projection, the error correction terms are 
re-scaled correspondingly down to their correct magnitudes 
by varying the exposure period. As the number of iterations 65 
increase, the magnitudes of the errors decrease. Therefore, 
the distortion caused by the noise associated with the optical 



transform and the finite dynamic range of the device are 
minimized by reducing relative errors using the feedback 
scheme. 

FIGS. 46a-d show the result obtained by simulating the 
IFBP algorithm for synthesizing a 128x128 image. The 
dynamic range of the SLM's and the filter are set at 256 and 
±16 respectively. The system is ovexdetermined (M>N) to 
the same degree as that employed for the simulation of the 
SIRT. The result shows that the first reconstructed image is 
distorted due to the poor accuracy associated with the 
filtering. The distortion is corrected after the second and 
third iteration. When these results are compared with the 
results obtained in the simulation of the SLM/CCD imple- 
mentation of SIRT in FIG. 24, it is seen that the number of 
iterations required to obtain the same level of performance 
is significantly lower when the IFBP method is employed. In 
addition, the distortion is reduced when the numbers of 
distinguishable dynamic ranges of the SLM's are set at the 
same level as before, i.e., 256. 

FIGS. 41a-d show reconstructed images (32x32) con- 
verging to the original cross-section. The number of projec- 
tions is 60 and projections are taken from a 32x32 original 
cross- sectional image. 

The IFBP algorithm is an iterative version of the FBP 
algorithm. With the IFBP, iterative reconstructions are 
speeded up significantly and distortion in the reconstructed 
images is reduced. When the system is overdetermined, the 
solution of IFBP minimizes the least square errors of the 
filtered projections in the frequency domain with a square 
root radius filter. 

The potential application of this algorithm is in high speed 
and low cost optoelectronic systems. Inaccuracies associated 
with optical systems is minimized within a few iterations. 
The structure can also be used to implement the Projection 
Space Iterative Reconstruction-Reprojection (PIRR) algo- 
rithm in limited view tomography applications. 

Accelerated Method Using a FIR Filter 

In this section, another accelerated method utilizing an 
FIR filter is proposed. In addition, a hybrid prototype using 
off-the-shelf devices for implementing the algorithms is 
proposed. 

The convolution back projection method (CBP) for direct 
tomographic image reconstruction utilizes spatial domain 
FIR filtering instead of frequency domain filtering for the 
filtered back projection method (FBP). An hybrid structure 
implementing digital FIR filtering and optical back projec- 
tion has been proposed by Gmitro, et. al. However, the large 
size of the FIR filter increases the cost and processing time. 
In addition, distortions caused by the optical back projection 
stage cannot be avoided. The PERT iterative procedure, not 
only allows distortion to be corrected, but also the length of 
the FIR filter can be significantly reduced. It will be shown 
that the reconstruction can be accelerated process consider- 
ably by using only a few taps. It will also be shown that the 
PIRT with partial convolution (PIRT-PC) converges to the 
minimum-norm solution for an underdetermined system. In 
order to circumvent the limitations of finite dynamic range 
of the both back projection and forward projection SLM's, 
an alternative iterative procedure — called the residual itera- 
tive reconstruction technique (RIRT) is presented. 

Hybrid Approach 

The rationale behind optical implementations of tomo- 
graphic forward and back projections has been described 
extensively in literature. The optical forward projection 
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processor is also called optical Radon Transformer and 
employs the same principle as the optical Hough transform 
feature space processor. Optical back projection processors 
have been used for implementing direct reconstruction algo- 
rithms. The optoelectronic implementations of iterative 
algorithms using CCD and SLM arrays were discussed 
previously. A hybrid prototype consisting of an SLM/CCD 
back projection processor, a digital auxiliary image 
processor, an SLM/CCD forward projection processor and 
an auxiliary projection processor is proposed in the section. 
The SLM and the CCD arrays in the back projection 
processors are a linear FLC SLM (Ferroelectric Liquid 
Crystal SLM) array and a commercial CCD image detecting 
array respectively. The SLM and CCD arrays in the forward 
processor consists of a LCTV SLM (Liquid Crystal TV 
SLM) panel and a linear CCD detector array respectively. 
The auxiliary processors in the image and projection spaces 
can be built using commercial microprocessors. The two 
additional electronic digital microprocessor-based auxiliary 
processors are important for flexibility and ensuring stabil- 20 
ity. 

In an analog electrically addressed SLM, the output light 
intensity is proportional to the applied voltage. Several 
distinguishable gray levels can be obtained using such 
SLM's. In the case of binary SLM's, the output can be either 
in the "ON" or "OFF" states. Analog SLM's are very often 
used in image and signal processing applications. However, 
commercial available analog SLM's, such as the LCTV, can 
only operate at the speed of about 100 Hz. Binary SLM's, 
such as FLC SLM, can however be operated at a rate of 30 
about 100 KHz. In many applications, analog SLM's oper- 
ating at KHz rates are desired. In this section, a method of 
implementing high speed analog SLM's for generating 
multiple gray levels using fast binary SLM's is described. 

The convolution back projection method can be derived 
from the filtered back projection method. A filtered projec- 
tion in eq. 75 can be described as 
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where H(p) = 



IpKpo 
otherwise 
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and P(0, r) is the Fourier representation of the measured 
projections. The frequency response is shown in FIG. 38. 
The impulse response, h(r) in the spatial domain, is given by 
the inverse Fourier transform of H(p) which is 



[sin27i/po i / sinrpo \ 1 
2Tt^o "2\— — ) J 



(92) 
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Using the convolution theorem, the filtered projection in 
eq. 91 can be written as 

f . < 95 > 

Q(6,r) = J *>(e,Ofc(r- 

where b(6, r) is the measured projection. In discrete form 

2(0,0= £ h&ybi&j-k) (96) 



As in the case of the FBP method described by eq. 76, the 
cross-sectional image can be reconstructed by superimpos- 
ing all filtered projections at different projection angles in 
the image plane. Using the convention defined in eq's. 91 to 
93, the convolution back projection method can be 
expressed in matrix form. We truncate h[i] and incorporate 
a finite number of samples in the impulse response the 
matrix H. Then, the reconstruction scheme can be expressed 
as 

x=A T Hb (97) 

where H is an MxM matrix representing the zero-padded 
convolution operation given by 
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and m is the number of ray-sums per projection. The matrix 
H is not only symmetric but also positive definite because it 
is diagonally dominant for any finite m. The diagonal 
dominant nature can be shown as follows, 
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= po 



A(0)+.£ lft(f)l = po 
r=-c« 



In the discrete case, we substitute the sampling interval T=V2 60 

Po in eq. 92, to obtain the impulse response h(i), i=-oo to oc, Therefore, when m is finite, 



r./-\ o F simri l / siruri/2 ^ 1 
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h(0)> Z ft(i)l 
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which is shown in FIG. 48. The impulse response can be It has been shown that, when the dimensions of the system 

rewritten as are finite, H is symmetric and positive definite and the 
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inverse, H _1 , exists. A reduced order convolution scheme and 
can be used as an FIR filter to accelerate the convergence of 

the PIRT algorithm. The improved algorithm can be x„- T, aA T Hn ^ 106 ^ 

expressed as ^ 



^jfc-i-HX H [b-A x k _,\ 



Substituting eq. 105 into eq. 106 



x k ^A T t k (101) *n = (X/tf?/{ ^[I-ctAATHf J «> 

where p=l/||AA 7 ]| 2 is used to normalize the projections and 10 = vA T H{i -[/- uAA T ff\}- 1 

a is chosen so that ||I~-apAA :rr l| 2 <l. The solution will f , n . ^ 

converge to the minimum-norm solution provided the sys- | {/-[/- txAA 7 ?/]} j - ctAA T ff] k J J m 

tern is underdetermined 1 
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When ||I~aAA :rr || 2 <l, and n->°o. 
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P The solution, x n , converges to 
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where r 0 =b 

25 An alternative iterative procedure can be written as 



Eq. 101 represents the basic PIRT accelerated by partial 
convolutions (PIRT-PC). The coefficients of the n th order r^r^-a A A T H 

filter are given by t k =t^+a H r k (110) 

p 0j i = 0 (103) 30 The final solution is 

-po - \ ~ » z odd and i ^ « (111) 
7t 2 i 2 . T 

* = A* Hm 

0, otherwise. 

Employing the reduced convolution, reconstructions of the 35 For an under determined system, the solution obtained using 
high frequency components are speeded up. Usually, only a tne algorithms described by eqs. 104 and 110 are the same 
few taps are needed in order to accelerate the convergence if numerical errors are ignored. If the system is 
of high frequency components. When the order of the FIR overdetermined, the state vector t given by eq. 100 will 
filter reduces to zero, H becomes an identity matrix reducing 40 diverge in the projection space, 
me approach to the basic PIRT. Hybrid Solution 

Let r*. denote the residual in the projection space. A 
tomographic reconstruction system using the information Optoelectronic architectures for implementing a variety 

contained in the residual to arrive at the minimum-norm of rt f ative reconstruction algorithms have been outlined in 
solution can be described as follows 45 previous sections This section describes a hybrid prototypic 

system that can be built to establish the proof of principle. 
aAA r /fr* i ( 104 ) ^ e hybrid system employs off-the-shelf devices. 

FIG. 49 shows a line diagram outlining the basic building 
r*_i + k-i blocks in the system The building blocks include 

dk = uA T Hn 50 An Optoelectronic Back Projection Processor. 

xk = x^x + dfr An Optoelectronic Forward Projection Processor. 

TT . . . t . , An Image Space Auxiliary Processor, 

where H is the partial convolution matrix and a is chosen so . . « ^ 

mat||I-aAA^H||<l.Theinitialvaluescanbechosenasr 0 =b £ A ^ h ^ ^T^u ^ « ■ * 

. rr-rf j j , . . ~ . . . M + " c ^ Each of the building blocks is described briefly m the 

and Xq=0. The method described by eq. 104 is similar to the 55 f Q n ow £ n sec tions 

steepest descent method with a fixed descent step a. The VhTbtck projection operations in tomographic image 
matrix H can be considered as a preconditioning matrix. The reconstruction involves mapping the data in the projection 
algorithm expressed by eq. 99 will be referred to as the space mto ^ imagG space from ^ ±c measured projection 
residual iterative reconstruction technique (RIRT). angles. The optoelectronic back projection processor shown 
The convergence properties can be derived as follows. By in FIG. 50 performs the tomographic back projection opera- 
substituting r k _ m by [I-ocAA 77 ] r£_ (m+1) into eq. 104 recur- tions. It consists of: 

sively until m=k-l, we obtain a linear array of analog SLM's and a cylindrical lens. 



= [I-oAA T H\n-\ ( 105 ) 



an image detection array. 
65 an image rotator. 

[/ - cfAA T H] k m The transmissivity of the linear SOvl array represents the 

magnitudes of the data in the projection space corresponding 
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to a projection angle. The pattern in the linear array is then 
"stretched" into an array of strips as shown in FIG. 50. The 
linear SLM array is not specified deliberately since a wide 
range of choices are available. Among those choices is the 
Ferroelectric Liquid Crystal (FLC) SLM.. 5 

The image detection array integrates the projection data, 
corresponding to the area intercepted on the pixel by the 
light beam, at all the projection angles. The image detection 
array can be a CCD image detector or a commercial CCD 
TV camera. 10 

A complete set of back projections involves performing 
back projection operations at all the projection angles over 
an angle of 180°. The purpose of the image rotator is to 
rotate the projection pattern to the corresponding projection 
angle relative to the image detector. The most commonly 15 
used image rotator consists of a dove prism image rotator 
driven by a stepper motor. FIGS. 51a^c show a dove prism 
image rotator where the prism is rotated 0, 45 and 90 degrees 
and the image is rotated 0, 90, and 180 degrees correspond- 
ingly. 20 

The dove prism image rotator is not electronically con- 
nected to the mechanical rotating part. However, optical 
distortions may be introduced by the prism. The operation 
can also be performed by rotating the linear SLM array and 
the lens "back and forth" over a range of 180°, Alternative 25 
approaches using a static structure can also be employed. 

The Optoelectronic Forward Projection Processor 

The optoelectronic forward projection processor as shown 
in FIG. 52 performs the forward projection operation to map 30 
the data in the image space into the projection space. It 
consists of: 

an image display array. 

a linear projection detector array and a cylindrical lens. 

an image rotator. 

The image display array displays a reconstructed image 
using either a two dimensional analog SLM or a commercial 
liquid crystal projection TV. The linear projection detector 
array detects the image generated by the cylindrical lens. 4Q 
The linear projection detector array can be a linear CCD 
array. The role and the design of the image rotator is 
identical to the one described in the preceding section. 

Since operations of image addition can be performed 
directly on the CCD array, a separate image auxiliary 45 
processor may not be required. However, when the perfor- 
mance of devices is not guaranteed, additional flexibility and 
improvement in accuracy can be gained by employing an 
auxiliary processor. The processor is required to be able to 
store a frame of the image and perform simple pixel by pixel 5Q 
addition or subtraction. Real-time digital video image pro- 
cessors which are commercially available can be used to 
perform these operations. 

The basic functions of the auxiliary space processor are: 

to receive data from the forward projection detector 55 
arrays. 

perform simple arithmetic operations, 
store projections, 
and 

send data to the back projection linear SLM array. 60 
This basic auxiliary processor can be built using either a 

commercial digital microprocessor or a special purpose 

processor. 

Most commercial analog SLM's are developed for build- 
ing television displays. Examples include liquid crystal 65 
SLM panels used for image projectors. The refresh or update 
periods for analog SLM's are typically in the range of video 
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frame rates, i.e., one to ten milliseconds. Unfortunately, the 
algorithms described earlier involve several hundreds of 
projections in each iteration. If such analog SLM's are used 
for implementing the 1-D back projections, then the result- 
ing processing period will be in the range of a few seconds. 
In order to increase the processing speed, high speed 1-D 
analog SLM's are required. This can be accomplished by 
controlling the duration of exposure instead of relying on the 
need for varying the gray levels continuously. 

CCD detectors are essentially integrating devices gener- 
ating charge/voltage levels proportional to the intensity of 
the incident light and the exposure period. When the CCD 
array is coupled with a SLM array, the charge/voltage levels 
in the CCD can be controlled by varying either the intensity 
of the light being integrated or the "ON" period of the SLM. 
Since most fast SLM's are binary in nature, we proposed to 
use the latter approach to control the CCD signal. The timing 
diagram of scheme is shown in FIGS. 53a-b. 

This circuit consists of a comparator, an address counter, 
a timing counter, timing memories, and registers as shown 
in FIG. 54. Values of grey levels of the SLM cells are 
preloaded into the timing memories. A cell is switched 
"OFF" when the value of the timing counter is larger than 
the value stored in the corresponding timing memory. The 
cell remains "ON" otherwise. During each timing counter 
clock period, the contents of the timing memories are read 
out sequentially with the address generated by the address 
counter. 

An example involving four SLM cells and sixteen gray 
levels is shown in FIG. 55. The address counter is clocked 
at a period T address- Th e timing counter is clocked by the 
carry output of the address counter whose period is r T timing . 
The period T^ n ^=mxT aAfms . s . where m is the number of 
SLM cells driven by the circuit. During each T ftmi „^, all the 
timing parameters in the memories are scanned out and 
compared with the output of the timing counter. If the value 
read out from the memory is less the value of the counter, 
then the corresponding SLM cell is set to one, else it is reset 
to zero. The maximum number of gray levels which can be 
represented is given by 



T- address X ttl 



(112) 



where T toto/ is the frame period for an analog SLM and 
^update is the time required for updating the memories. When 
the number of SLM cells is large, the speed depends on the 
speed of the circuit at which it can be driven. 

After all the n gray levels are counted, the memories are 
updated by the external data with a fast memory access 
technique, such as a direct memory access scheme (DMA). 
When a large number of SLM cells are to be driven, the cells 
can be partitioned into blocks. The size of blocks depends on 
the desired driving speed. 

When the timing clock period is less than the SLM 
settling time, the resolution deteriorates. The problem can be 
circumvented by setting a constant zero bias exposure time 
which is equal to the minimum SLM settling time. The 
additional offset can be subtracted from the output of each 
CCD detector. 

The dynamic range of reconstructed image is mainly 
affected by three factors: 

a) The nonnegative offset associated with the projection data 
may reduce the actual dynamic range to approximately 
10%. 

b) The finite dynamic range of the back projection 1-D SLM 
array. 

c) The finite dynamic range of the forward projection 2-D 
SLM array. 
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The first effect can be eliminated by back projecting the 
positive and negative projection data separately. In addition, 
the error caused by the finite contrast ratio can also be 
rninimized by using the dual projection scheme. 

The RTRT described in eq. 104 provides a method for 
minimizing distortions caused by the finite dynamic ranges 
of the SLM arrays. In the iterative procedure, only back and 
forward projections of residuals are involved. In a stable 
feedback system, the magnitude of the residual will vanish 
as the number of iterations increase. Therefore, the vanish- 
ing residuals can always be adjusted to fit the maximum 
dynamic range of the SLM arrays and then rescaled back 
appropriately using the auxiliary processors after optical 
projections. 

FIGS. S6a-d show the reconstructed Shepp and Logan 
phantom at the 8"*, 16 th , 32 nd . and 64 th iterations when a 
third order FIR filter (n-3, n-1, n, n+1, n+3) is used. FIGS. 
S7a—d show the reconstructed images obtained using the 
basic PIRT at the same iterations. FIGS. 58 and 59 are 
cross-sections of reconstructed images using PIRT- PC and 
PERT respectively where the number of iterations are 4, 8, 16 
and 32. It is apparent that the convergence of high frequency 
components is accelerated considerably when a short FIR 
filter is used. However, the speed of convergence of the low 
frequency components does not improve. 

FIGS. 60 and 61 show plots of the cross sections at the 
32 nd and 64 th iterations respectively. The error across the 
image, as shown in FIG. 60 is less than 1% of the entire 
dynamic range except for regions where the variation is 
sharp. The plot shown in FIG. 61 indicates that there is very 
little improvement in the result after 32 iterations. If com- 
mercial video devices operating at 30 frames per second are 
employed, then the reconstruction time, with an expected 
contrast resolution of less than 1%, is about 1 second 

Finite dynamic ranges (numbers of distinguishable grey 
levels) of both forward and back projection SLM arrays 
represent the dominant source of distortion in the structure. 
The dynamic ranges of the available CCD arrays are much 
higher than those of available SLM arrays and consequently 
the effects can be neglected. In the simulation, 256 distin- 
guishable gray levels are imposed by normalizing the state 
variable t and the reconstructed image x and quantizing them 
into 255 distinguishable levels before projecting. In 
addition, the reconstructed image is constrained to be within 
the range 0 to 255. FIGS. 62a-d show reconstructed images 
at the 8'*, 16'*, 32 nd and the 64'* iterations. The results show 
that, compared to the image obtained by PIRT- PC without 
quantization as seen in FIGS. 56a-<i, the distortion is 
insignificant. 

If all projections are biased to nonnegative values, then, 
only about ten percent of the dynamic range in the image 
space is available. FIG. 63 shows a plot of the cross- section 
of a reconstructed image using this strategy. However, if 
projections are separated into positive and negative projec- 
tions and back projected separately, then, there is no deg- 
radation in the useful dynamic range in the image space as 
shown by the plot in FIG. 64. 

The PIKT-PC proposed in this section is a combination of 
PIRT and reduced CBP. This algorithm is capable of reduc- 
ing the cost associated with a full size convolution with CBP, 
correct the errors of reduced length convolutions, and pos- 
sesses superior convergence properties relating to PIRT. In 
addition, the method guarantees the minimum-norm solution 
for an underdetermined system. 

The proposed hybrid prototype described in this section 
represents a practical design which can be built using 
off-the-shelf video and optoelectronic devices. The proto- 
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type design represents a first step. The additional electronic 
processors in the projection and the image space not only 
offer flexibility which is necessary while building the first 
prototype. But also enables us to investigate other algo- 
rithms. In comparison with existing approaches, the pro- 
posed structure offers 

* reduced dependency on the accuracy of filter devices. 

* reduced dependency on the precision and dynamic range 
of optoelectronic devices. 

* significant reduction in the cost of filtering devices. 
Several iterative tomographic image reconstruction meth- 
ods have been presented. These methods are based on the 
projection iterative reconstruction technique. In contrast to 
conventional approaches, the algorithms described herein 
offer superior results when the tomographic system is under- 
determined and computational errors are present. 

New system concepts, involving the use of optoelectronic 
devices for implementing iterative tomographic image 
reconstruction algorithms, have also been proposed. The 
optoelectronic feedback structures built with conventional 
optical and imaging devices offer high speed reconstructions 
with quality levels required in medical applications. 

The PIRT algorithm has been developed for tomographic 
image reconstruction. As shown in Table 5, the new algo- 
rithm differs from the conventional iterative algorithms, 
such as SIRT, in that it solves for a state vector in the 
projection space in contrast to conventional algorithms 
which solve for the solution in the image space. The feed 
back process in PIRT is performed in the projection space. 
Hence, iterative error corrections in the image space are 
avoided. The feedback scheme in SIRT is implemented in 
the image space leading to complications in optoelectronic 
implementations. 

TABLE 5 

SIRT and PIRT Algorithms 



algorithm SIRT 
expression x(k) = x(k — 1) + A T (b - 

Ax(k -1 )) 
feedback image space 



PIRT 

t(k) = t(k - 1) + (b - AA T t(k - 
1)) 

projection space 



When a system is underdetermined, the new algorithm PIRT 
is superior to the conventional SIRT algorithm since it 
guarantees a unique minimum-norm solution. The properties 
of the two algorithms are compared in table 6. 

TABLE 6 

Comparison of Performance for Underdetermined Systems 



algorithm 


SIRT 


PIRT 


solution 


not unique 


minimum-norm 


convergence 


not unique 


guaranteed 


eg method 


does not converge 


guaranteed 


initial 


no unique solution 


minirrmm norm of 


image 




difference 


with 


no unique solution 


minimum-norm in reduce 


constrant 




space 



65 However, the PIRT is not useful for solving overdeter- 
mined systems since the state vector may diverge in the 
projection space as show in table 7. 
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TABLE 7 



Comparison of Performance for Overdetermmcd Systems 

algorith SIRT PIKT 

solution least squares least squares 

convergence guaranteed diverges in projection 

space 



The tomographic reconstruction algorithms described in 10 
this dissertation include PIRT, PIKT-CG. PERT-PC, RIRT 
and IFBP. 

The PIKT is a basic iterative algorithm for solving large 
sparse underdetermined systems. The PIRT is considered as 
a counterpart of the conventional iterative algorithm SIKH 15 
The PIKT solves for the minimum-norm solution of normal 
equations whereas the conventional algorithms are based on 
a strategy for computing the least squares solutions. The 
PIRT guarantees unique solutions for underdetermined sys- 
tems while conventional algorithms provide unique solu- 
tions when systems are overdetermined. 20 

The PIRT-CG method differs from the conventional SERT- 
CG methods in that it applies the conjugate gradient accel- 
eration method to solve for a state vector in the projection 
space instead of directly solving for a solution vector in the 
image space. Therefore, when the system is 25 
underdetermined, the algorithm does not diverge since the 
matrix of the state equation is not singular. 

The PIRT-PC method is another accelerated PIRT type 
algorithm in which a short order FIR filter is used to speed 
up the convergence of high frequency components. Results 30 
show that when a third order filter is used, the convergence 
can be accelerated significantly. The PIKT-PC is a simple 
and efficient algorithm, and should therefore be preferred for 
optoelectronic implementations. 

The RIRT recursively eliminates the residue in the pro- 35 
jection space instead of building a state vector. RIRT algo- 
rithms are equivalent to SIRT, PIRT, PIRT-CG, PIRT-PC or 
the steepest descent type algorithms depending on the accel- 
eration technique used. When computational errors exist, the 
algorithm does not guarantee uniqueness of the solution. 40 

The EFBP is a SIRT type algorithm whose convergence 
properties are improved by using a spatial frequency filter. 
This algorithm is mainly useful for optoelectronic imple- 
mentations since optical devices can be used to build a fast 
spatial frequency domain filter although the accuracy is 45 
relatively poor. The structure uses a feedback scheme to 
eliminate distortions associated with optical transforms. 

The structure implementing SIRT is simple and does not 
employ sophisticated filtering or convolution processing 
units. The structure for implementing IFBP using optoelec- 50 
tronic devices employs a high speed, low accuracy, inex- 
pensive optical filter to accelerate convergence. 

In the case of SIRT, the error corrections are superim- 
posed on the previously reconstructed image. Therefore, 
errors introduced by optoelectronic devices may be accu- 55 
mulated in the reconstructed image. This is especially true, 
if the system is underdetermined where the accumulated 
errors may cause divergence and reach the dynamic range 
limits of devices. With the proposed PERT algorithm, the 
need for superimposing images is eliminated. This leads to 60 
a smoother image with a better spatial resolution from fewer 
projection data. The PERT has reduced the complexity of 
structures for optoelectronic tomographic reconstruction. In 
addition, schemes for acceleration, such as PIRT-PC, can be 
implemented without additional expense. 65 

Iterative tomographic image reconstruction algorithms 
are commonly used in single photon emission computerized 
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tomography (SPECT) and positron emission tomography 
(PET). Applications of PIRT type algorithms for SPECT and 
PET applications need to be further developed. 

Many iterative estimation algorithms have been devel- 
oped based on statistical models of tomographic systems. 
Some of these algorithms can be implemented using PIRT as 
a tool. 

Only optoelectronic implementations of additive iterative 
reconstruction algcrithms have been presented. Multiplica- 
tive iterative reconstruction algorithms are also commonly 
used. Multiplicative algorithms offer solutions which are 
nonnegative and are therefore attractive for optoelectronic 
implementation. 

Multiplicative x-ray tomographic algorithms can be 
implemented with a slightly modified optoelectronic struc- 
ture. In the modified structure, only an additional SLM array 
needs to be coupled to the back projection CCD array. The 
system concept is straightforward. However, the theoretical 
foundations are yet to be developed. 

A hybrid prototype for implementing the iterative algo- 
rithms has been proposed. The hybrid structure can provide 
flexibility required for pursuing further research. The pro- 
totype can be built using commercial available optoelec- 
tronic devices and reach a speed of 30 frames per second. 

Novel optoelectronic systems for implementing iterative 
tomographic image reconstruction algorithms have been 
proposed. In order to exploit the advantages of optoelec- 
tronic devices fully, a new iterative tomographic image 
reconstruction algorithm — PIRT has been developed. The 
PIRT algorithm is general, in that it can be employed to 
solve a wide variety of problems involving large sparse 
matrices and underdetermined systems. The ideas developed 
in the present invention are, therefore, far broader in scope 
than their narrow application to tomographic reconstruction 
described here would suggest. 

Summary 

The tomographic reconstruction scheme outlined in this 
disclosure offers many advantages over conventional meth- 
ods of implementation. In comparison with existing 
approaches, the method outlined herein offers 

reduced dependency on the accuracy of filter devices. 

reduced dependency on the precision and dynamic range 
of optoelectronic devices. 

significant reduction in the cost of filtering devices. 

The new algorithm PIRT is superior to the conventional 
algorithm SIRT since a unique solution is guaranteed when 
the system is underdetermined. The PERT-PC is a combina- 
tion of the PIRT and reduced CBP. This algorithm is able to 
reduce the cost of a full size CBP, correct the errors of 
reduced CBP, and accelerate the convergence of the PIRT. 

We claim: 

1. A back projection processor, comprising: 

(a) a 2-D charge coupled device (CCD) detecting array; 

(b) spatial light modulator (SLM) projection means for 
generating, using spatial light modulation, in response 
to input data representing a 1-D image, a corresponding 
1-D image stretched along one axis to provide a 
stretched 1-D image, and for projecting the stretched 

1- D image onto the 2-D CCD detecting array; 

(c) image rotation means for orienting a stretched 1-D 
image projected by the SLM projection means on the 

2- D CCD detecting array at a desired projection angle, 
the 2-D CCD detecting array outputting back projected 
data; and 

(d) whereby the back projection processor smears a 1-D 
image represented in a data array onto the 2-D CCD 
detecting array. 
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2. A forward projection processor, comprising: 

(a) a CCD detecting array; 

(b) SLM projection means for generating, using spatial 
light modulation, in response to input data representing 
a 2-D image, a corresponding image compressed along 
one axis to provide a compressed 1-D image, and for 
projecting the compressed 1-D image on the CCD 
detecting array, the CCD detecting array outputting 
forward projected data, 
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(c) image rotation means for orienting a compressed 1-D 
image projected by the SLM projection means on the 
CCD detecting array at a desired projection angle; and 

(d) whereby the forward projection processor forward 
projects a 2-D image represented in a data array onto 
the CCD detecting array. 
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